Three Dimensional Numerical Relativity with a Hyperbolic Formulation 
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We discuss a successful three-dimensional cartesian implementation of the Bona-Masso hyperbolic 
formulation of the 3+1 Einstein evolution equations in numerical relativity. The numerical code, 
which we call "Cactus," provides a general framework for 3D numerical relativity, and can in- 
clude various formulations of the evolution equations, initial data sets, and analysis modules. We 
show important code tests, including dynamically sliced flat space, wave spacetimes, and black hole 
spacetimes. We discuss the numerical convergence of each spacetime, and also compare results with 
previously tested codes based on other formalisms, including the traditional ADM formalism. This 
is the first time that a hyperbolic reformulation of Einstein's equations has been shown appropriate 
for three-dimensional numerical relativity in a wide variety of spacetimes. 



I. INTRODUCTION AND OVERVIEW 

The young field of three-dimensional (3D) numerical 
relativity has entered an exciting era. As we review be- 
low at some length, strong theoretical and astrophysical 
motivations have led to increased activity and collabora- 
tions among many research groups, and general 3D rela- 
tivistic problems are being attacked with increasing suc- 
cess. In this paper we present results from a new and 
very general 3D advanced computer code, which we call 
"Cactus" , designed to study these general problems in 
a collaborative environment. This is the first in a series 
of papers on this code and its applications in numeri- 
cal relativity and relativistic astrophysics. At the same 
time, the paper is also a follow-up of previous papers in 
our continuing exploration of hyperbolic formulations of 
the Einstein equations for numerical relativity B-ff] . 



A. Motivation 

The imminent arrival of data from of the long awaited 
gravitational wave detectors (LIGO, VIRGO, GEO600, 
TAMA; see, e.g., Rcf. §| and references therein) has pro- 
vided a sense of urgency in producing realistic simula- 
tions of very strong sources of gravitational waves, which 
can only be done through the full machinery of numerical 
relativity. One of the best candidates for early detection 
by the laser interferometer network is increasingly consid- 
ered to be black hole mergers |||5). However, the signals 
are likely to be weak enough by the time they reach the 
detectors that reliable detection may be difficult with- 
out prior knowledge of the merger waveform. These are 
among the reasons that the NSF-funded Binary Black 
Hole Grand Challenge Alliance has focused the efforts 
of numerous US and international groups on developing 



codes for solving the problem of 3D coalescing black holes 
(see, e.g, the latest round of papers of the Alliance ||-§). 

Another important process in astrophysics that re- 
quires fully relativistic simulation is neutron star mergers 
(see, e.; 
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which will produce a possibly detectable 
burst of gravitational waves [ pd| . These are sometimes 
considered as sources of gamma-ray bursts [[l2|, and the 
final state (e.g., a neutron star or black hole) is highly 
uncertain. Most studies of this process have been New- 
tonian, but even post-Newtonian correction terms, which 
are still inadequate to describe the possible formation of 
a black hole, produce significant changes in the evolu- 
tion |l3] ]. More relativistic approximations to the Ein- 
stein equations produce still quite different (and contro- 
versial) results, indicating that the neutron stars may 
actually form black holes before the merger Jul. The 
point we wish to make is that the merger process clearly 
requires a fully consistent relativistic treatment, which 
provides another motivation for development of powerful 
and general numerical codes to solve the full set of Ein- 
stein equations, in this case coupled to the relativistic 
fluid equations. This research area is a particular appli- 
cation for which Cactus is being developed, although we 
will present only tests of the vacuum part of the evolution 
system in this paper. 

Astrophysics aside, there are of course purely theoret- 
ical reasons to develop robust 3D solvers to Einstein's 
equations. As general relativity is one the fundamental 
theories of physics, it needs to be better understood in its 
most nonlinear regimes, which are usually the most diffi- 
cult to probe. Again, numerical treatment of the full set 
of Einstein equations is one of the main tools for study- 
ing the theory in such regimes, and has already led to 
the discovery of unexpected phenomena, such as critical 
phenomena in black hole formation |14| (see Ref. p| for 
a recent review), which has now been seen in spacetimes 
containing scalar fields, fluids, and even in pure vacuum, 
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gravitational wave spacetimes. Most of these studies have 
been carried out in ID or in rare cases in 2D |f6f| , but 
little is known about the 3D behavior (l?]]. 

Unfortunately, despite all of these motivations for 3D 
numerical relativity, and the best efforts of many groups 
around the world, progress has been slower than hoped 
and expected. One of the reasons for this is the the sheer 
complexity of the Einstein equations in 3D, coupled with 
the immense computational needs for solving them. For 
example, an enormous amount of memory and time on 
the order of one CPU day on a Teraflop computer will be 
required to produce a single, highly resolved simulation 
of 3D black hole spiraling coalescence (see Ref. for 
a review). Developing well tested software that simulta- 
neously solves the Einstein equations, takes advantage of 
high performance parallel computers, and can be effec- 
tively used by the large number collaborators needed to 
develop algorithms is a challenging software engineering 
problem in its own right. 

However, the problems of 3D numerical relativity run 
far deeper than computation and code development. 
Given a sufficiently large computer and perfectly de- 
bugged code, problems like coalescing black holes or neu- 
tron stars would still not be solvable today, because of 
important theoretical and algorithmic problems still to 
be addressed. Perhaps the best example to illustrate 
these problems is that of a spacetime containing black 
holes. 

The presence of a singularity inside the black hole and 
the weak field zone far from the hole gives rise to an 
extreme dynamic range. Singularity avoiding slicings ef- 
fectively keep time slices from hitting the singularity, but 
lead to pathological time slices that create huge gradients 
near the black hole horizon which cannot be resolved, es- 
pecially in 3D [|l9|,^0|. Such gradients lead to numerical 
instabilities with the standard formulations of the equa- 
tions, often causing codes to crash in the interior well 
before the desired evolution can be carried out in the ra- 
diation zone. Although the characteristic time necessary 
to obtain accurate waveforms for the inspiral and merger 
of two black holes is on the order of thousands of M (see, 
e.g, Ref. J9|), even state-of-the-art black hole collisions in 
axisymmetry (2D) [ pl|j22| can only be evolved for hun- 
dreds of M. (We will use units c = G = 1 throughout 
this paper, so time and spatial units for black hole sim- 
ulations are in terms of the black hole mass M.) 

Success in evolving black holes in 3D has been mixed. 
Partial successes include colliding, equal mass black holes 
in 3D [L8|. and waveform extraction of distorted 3D black 
holes [23 25 1 . In both cases the system is evolved suc- 



cessfully for tens of M and although this would be com- 
pletely insufficient for the black hole coalescing problem, 
it is enough time to study the waveforms produced in the 
ring-down phase. As verified by comparison with pertur- 
bation theory and axisymmetric simulations, these 3D 
simulations can produce highly accurate waveforms, but 
they ultimately crash both due to the "grid stretching" 
pq| effects created by singularity avoiding slicings and 



due to poor outer boundary conditions. 

One approach to understanding the expected wave- 
forms that avoids these problems is to solve the linearized 
equations describing black hole gravitational wave inter- 
actions. This approach has proven to be remarkably 
robust in comparisons with a range of presently feasi- 
ble fully nonlinear simulations of distorted and colliding 
black hole spacetimes ^5| , p7| -^T|, but it does not solve 
the general coalescence problem. Related studies of a di- 
rect 3D integration of the perturbation equations show 
that even such a simple linear problem is very demand- 
ing, having inner and outer boundary difficulties [ |32"[ that 
can be overcome through the machinery of adaptive mesh 
refinement p3|-p5[. 

The problem of dealing with singularities, grid stretch- 
ing, and inner boundaries may be ultimately solved by 
the so-called AHBC (apparent horizon boundary condi- 
tions) ]3^-^0|,pt, which are basically ingoing conditions 
on appropriate quantities evolving near the black hole 
horizon coupled with appropriate gauge conditions. But 
other gauge problems may still lead to large gradients as 
coordinates are sheared and squashed during the evolu- 
tion. Hence much research into appropriate gauge con- 
ditions for such dynamic spacetimes is needed. Even in 
very weak wave spacetimes, gauge problems can cause 
numerical codes to develop pathologies and crash as co- 
ordinates evolve out of control pl|-^3]|. Recent devel- 
opments shed new light into the mathematical under- 
standing of these coordinate problems and gauge patholo- 
gies in general ijUpq] , Furthermore, in order to resolve 
the inner, strong field region near the black holes, the 
outer boundary is generally placed uncomfortably close 
to the hole, where spurious signals or reflections which 
propagate inward may be generated due to inappropri- 
ate boundary conditions, masking the true physics taking 
place in the interior (for a recent discussion, see Ref. [Q). 

Despite these difficulties, there has been considerable 
progress in evolving dynamical black hole spacetimes in 
the last year. Brugmann |ff6|| recently demonstrated that 
it is possible to see some form of gravitational radia- 
tion from numerically constructed true 3D black holes 
with spin and momentum. Unfortunately, these feasibil- 
ity studies seem to indicate that current techniques have 
face more severe difficulties with these highly dynamical 
systems, and cannot yet provide useful information for 
realistic gravitational wave astronomy The Grand 

Challenge Alliance has developed outer boundary con- 
ditions which appear to allow accurate outgoing wave 
boundary conditions in three dimensional numerical rel- 
ativity 0. Moreover, using causal differencing and a 
careful inner boundary treatment, the Alliance has been 
able to transport a black hole several black hole radii 
across a grid || . Work by Daues and collaborators has al- 
lowed single black holes to be evolved beyond f 00M using 
dynamically determined gauges fl40| . However, none of 
these treatments have shown, to date, the ability to pro- 
duce a long time stable evolution for colliding or highly 
distorted black holes in three dimensions, and many dif- 
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ficult problems remain to be solved. 

Another very recent approach to 3D black hole evolu- 
tion that completely avoids the problems of grid stretch- 
ing is characteristic evolution, which has successfully 
evolved 3D rotating and distorted black holes for es- 
sentially unlimited time periods (t « 60,000M 
These spectacular results are achieved by using an ingo- 
ing characteristic foliation of the black hole spacetime, 
using the horizon as an inner boundary. However, it is 
not clear yet if this method will be viable for evolution 
of very highly distorted or colliding black holes, where 
focusing of ingoing light rays may create caustics, lead- 
ing to a breakdown of the foliation. Also, ironically, the 
method is presently most successful when a black hole is 
present, creating an S 2 x R topology; dealing with the 
so-called r = problem is difficult for any formulation of 
the Einstein Equations, and is avoided by using cartesian 
grids in the standard 3+1 formulations, but the charac- 
teristic method cannot use cartesian grids, and would 
therefore have to face this problem in the absence of a 
black hole (e.g., for the coalescence of neutron stars). 
Nonetheless, the possibility of very long time evolutions 
demonstrated with the characteristic evolution scheme is 
an exceptionally significant achievement that seems likely 
to provide an alternate and superior approach for an in- 
teresting class of 3D black hole spacetimes. 

B. Hyperbolic Numerical Relativity 

In recent years, much renewed research into theoret- 
ical foundations of numerical relativity has led to the 
development of hyperbolic formulations of the Einstein 
equations for numerical relativity, which have numerous 
advantages over the standard ADM formulation [Q . We 
have addressed in detail this issue in a previous publica- 
tion in this series [Q]. In summary, they (a) provide a 
much better starting point for the mathematical analysis 
of well-posedness and existence of solutions fi(|[50|| , (6) 
are better suited than the standard ADM formulation to 
modern numerical methods developed for computational 
fluid dynamics |5lJ] and promise to handle large gradients 
(c) are more adapted to providing natural bound- 
ary conditions either on the black hole horizon or at the 
outer edge of the simulation, and(rf) still allow a very gen- 
eral class of gauge conditions (many of which are yet to 
be developed) that will be needed to control coordinate 
motion (although see Ref. |j45) for caveats of hyperbolic 
choices in the gauge conditions). 

Reula has recently reviewed, from the mathematical 
point of view, most of the recent hyperbolic formulations 
of the Einstein equations |3(j) (This article, in the online 
journal "Living Reviews in Relativity" , will be periodi- 
cally updated). It is important to realize that the math- 
ematical relativity field has been interested in hyperbolic 
formulations of the Einstein equations for many years and 
some systems that could have been suitable for numeri- 



cal relativity were already published in the 1980's f32| , |53"|] . 
However, these developments were not recognized by the 
numerical relativity community until recently. 

Choquet-Bruhat and Ruggcri already commented in 
1983 [^2) on the possible importance of stable hyperbolic 
systems for numerical applications. Following this sug- 
gestion, Bona and Masso studied the numerical relativity 
implications of the harmonic slicing condition ]54j and 
the advantages of systems of balance laws from the nu- 
merical point of view [p5| . In 1992 they proceeded to de- 
velop the first hyperbolic formulation of the 3D Einstein 
equations with numerical relativity in mind Special 
emphasis was put on the idea of borrowing from the huge 
arsenal of numerical methods available from the compu- 
tational fluid dynamics community. 

A complete 3D code was developed with this formula- 
tion [^,^7), leading to an advanced parallel version de- 
veloped at NCSA called the "H" code. Different varia- 
tions on this code were used in numerous applications in 
relativity where it was extensively tested on pure wave 
spacetimes [ flT|] , and in computational science (see, e.g., 
p8| , p9| ). This code forms the basis for some of the tests 
presented here, and furthermore the computational sci- 
ence experience gained from developing this code was 
essential in developing the more powerful Cactus code, 
described below. However, this formulation was hyper- 
bolic only for harmonic slicing (which amounts to a sim- 
ple algebraic condition on the lapse: a oc ^J~g, where g 
is the determinant of the three-metric gij), and it did 
not consider a shift, making it suitable only for a limited 
range of problems in numerical relativity. 

For these reasons, the system was generalized to apply 
to an arbitrary shift and to an infinite family of lapse 
conditions, including maximal slicing, in which case a 
mixed hyperbolic-elliptic system results |KM. This sys- 
tem, currently known as the "Bona-Masso formulation" 
(BM), takes the flux conservative form, which already 
allows a wide class of modern numerical methods not 
possible with the standard ADM formulation, for any 
choice of lapse and shift. But it has the additional ad- 
vantage of being hyperbolic (i.e., diagonalizable) if the 
lapse is chosen from the particular infinite class of slic- 
ings defined below. This formulation showed its superi- 
ority over the standard formulation in spherical symme- 
try (ID) by evolving a black hole essentially indefinitely 
without apparent horizon boundary conditions. Due to 
the use of the eigenfields, the advanced numerical meth- 
ods available to such a formulation, and the improved 
outer boundary treatment afforded by the formulation, it 
was able to handle the large gradients that develop near 
a black hole with a singularity avoiding slicing. Details 
of these numerical techniques and boundary treatments 
are given in an accompanying paper in this series J6l| . 
We are presently working to carry these techniques into 
3D, and this paper takes the first step in addressing these 
issues. 

The BM system is now one among many hyperbolic 
systems, as other independent hyperbolic formulations of 
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Einstein's equations were developed J62|~|67| at about the 
same time as Ref. ||60| , To our knowledge, among these 
other formulations only the one originally devised in Ref. 
has been applied to spacetimes containing black holes 
, although still only in the spherically symmetry ID 
case (a 3D version using full AHBC is under development 

There is an additional important motivation for hyper- 
bolic systems in general relativity provided by the inter- 
est in relativistic hydrodynamics, which will be needed 
to study systems like colliding neutron stars. Traditional 
approaches to relativistic hydrodynamics treat the left 
and right hand sides of Einstein equations separately, 
with different numerical methods, independent update 
routines, and so forth. However, relativistic hydrody- 
namics has a single set of equations, mathematically and 
philosophically. If the entire set of Einstein equations, in- 
cluding the fluid equations (which should be considered 
as a subset of the Einstein equations) could be formu- 
lated as a single hyperbolic system, a unified numerical 
treatment of the entire system would be possible. 



C. Goals of this Paper 

For all of these reasons, it is essential to develop robust 
and general 3D numerical codes to attack the many prob- 
lems in general relativity and astrophysics waiting to be 
solved, testing and comparing the different formulations 
of the Einstein equations. With these strong motivations, 
this paper has a two- fold purpose: 

First, as follow-up of our previous work on the theo- 
retical basis of our formulation we present the first 
detailed testing of a hyperbolic formulation of Einstein's 
equations in 3D on a variety of spacetimes that have be- 
come established benchmarks for numerical relativity, in- 
cluding black hole and gravitational wave spacetimes. In 
this paper we will not try to advance the results of previ- 
ous 3D codes but we show for the first time that with 
standard numerical methods for balance law systems 
(MacCormack and Lax-Wendroff schemes, discussed be- 
low), the BM formalism compares well with the tradi- 
tional ADM formulation. In this paper we present results 
on the formulation in its most general form, allowing ar- 
bitrary slicings and shifts. This form does not allow for 
advanced numerical methods based on the eigenfields of 
a hyperbolic system, or advanced boundary treatments. 
Such methods are subject to further research and work 
is in progress to apply them to this system of equations. 
We also report on how to establish a set of techniques for 
rigorous verification and self-convergence testing. 

Second, we present a code, called "Cactus" , that pro- 
vides a general, high performance framework for 3D nu- 
merical relativity in a collaborative environment, allow- 
ing for a number of formulations of the equations, general 
gauge and initial conditions, different numerical meth- 
ods, analysis tools, etc. This code is being developed as 



a general tool to be used for many different problems in 
3D numerical relativity, such as those described above. 
The philosophy behind this approach is described in an 
accompanying paper 0. The performance and paral- 
lclization aspects are described in accompanying papers 
|fl^[rj| . Other tests of the code, including matter tests, 
horizon finders, waveform extraction, etc. will be pub- 
lished in future papers in this series, as a growing number 
of international collaborators are extending the capabil- 
ities of the current version. 

We proceed as follows: In Sec. [n] we discuss basic con- 
cepts of our code, including the systems of equations, 
coordinate systems, gauge choices, and numerical meth- 
ods. In Sec. Ill we discuss numerical issues, including 
methods, boundary conditions, and convergence testing. 
In Sec. |^ we treat dynamically sliced flat space mod- 
els to demonstrate simple yet powerful code tests. In 
Sec. |y| we focus on a series of weak gravitational wave 
spacetimes, replicating results from Ref. 0] and extend- 



ing their study to non-axisymmetric cases. In Sec. VI we 
treat black hole spacetimes with a wide variety of slicings, 
and compare with the analytical solution in the case of 
geodesically sliced black hole. In all our test cases, we 
obtain rigorous self-consistent convergence and, in those 
cases published before, excellent agreement with known 
results. 



II. THEORETICAL CONCEPTS 

In this section, we discuss some basic theoretical con- 
cepts and introduce the choices that we have imple- 
mented. We follow closely the BM theoretical formula- 
tion of Ref. Q . Some aspects of the ADM formulations 
are also discussed in others papers p9|,Bl|,pO|] . 



A. The BM formulation 

The BM formulation of the Einstein equations is dis- 
cussed in detail in a previous paper in this series [|J . For 
completeness, here we write the basic equations, although 
the reader is directed to Ref. Q for further details and 
discussions. One of the fundamental advantages of this 
formulation is that the whole system can be written in 
first order balance law form: 



d t vi + d k F>lu 



(1) 



where the vector u displays the set of variables, and both 
fluxes F h and sources S are vector valued functions. We 
stress that the fluxes F h and the sources S do not contain 
any derivative of the set of variables, which is crucial for 
analyzing the causal structure of the system and for the 
application of appropriate numerical methods. 
The vector u has the following 37 quantities: 



(gij , a, Kij , D kij , A k , V k ) 



(2) 
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where g^j, a, and Kij have their standard definitions. As 
we have introduced a first order system, the following re- 
lations act as algebraic constraints imposed on the initial 
slice only: 



(Bij + Bji)/a, 



(17) 



A k = d k lna, D kij = - d k c 



and the special combination 



Vi = Di 



D r 



(3) 



(4) 



is considered as an algebraic constraint which will hold 
if and only if the momentum constraint is satisfied |Q. 
We define D k j = g km D k ij, i.e., we use the three-metric 
gij to raise and lower indices on objects, even if they do 
not transform as tensors. This is just a notational conve- 
nience. We also note that the shift vector (3 1 is not in this 
dynamical set, as it is considered a given arbitrary func- 
tion whose spatial derivatives B k % = i d k f3 l are known at 
any time. 

The fluxes in the set of Eqs. (fy) are: 



F_ 



F k a 



0, 
0, 



FtKij = -I3 k + a [ D% - n/2 V k 9ij 
+1/2** (Aj+2Vj-D^) 
DJ) } , 



F_Dkij — 
FlA k 
F k Vi 



+ 1/2 5 k [A l + 2V l 
-(3 r D ri] + a (Kij - 
-p r A r + a Q , 
-(3 k Vi + B k - B k 



i) 



The sources for these equations are: 



S^ 9ij = -2 a [K 

S- a = 
S-Ka = 



2/3 r D 



-of Q + a(3 r A r , 



(5) 
(6) 
(7) 



(8) 
(9) 

(10) 



(11) 
(12) 



2{K ir BT + K ]r B? - KijB r r ) 
+a [ - (4) % - 2K k K kj + trKKi 



ri kj 



2D lk r D r k 



2D/ k D r k 



T k T 



kr ij 

(13) 



~(2D k ~A r )(D^ + D 3 l) 
+MVj - 1/2 Dj k ) + Aj{Vi - 1/2 D lk k ) 
+A 3 (V l -l/2 D k )-nV k D H j] 
+n/A a 9ij [ -D k rs T k rs + D k lD k \ - 2 V k A k 
+K rs K rs - (tr Kf + 2a 2 G 00 } , 
S-D ki j = , (14) 
S-A k = , (15) 
S^Vi = a [a G" + A r (K\ -trK 

+K r s (D rr s - 2D r !) - K\{D r s s - 2D J)] (16) 
+2{Bl - ^ tr B) V r + 2(D r ? - Sf D j jr )B r s . 

We have used the shorthand 



and we stress again that for notational convenience, we 
raise and lower indices with the three-metric, so for in- 
stance we have written B^ = g ik B k , even though Bij is 
not a tensor quantity. 

The free parameter n allows one to select a specific 
evolution system (it is zero for the "Ricci" system and 
one for the "Einstein" system) , as discussed in Ref. . 

As in this paper we do not explore methods based on 
the the diagonalization of the system (i.e., based on the 
characteristic fields), we will not detail here the spectral 
decomposition. The reader is directed to Ref. j| for all 
the theoretical foundations of hyperbolicity. Applications 
of advanced hyperbolic methods to the eigensystem in 
one-dimensional problems can be found in Ref. |6l| . 



B. ADM Formulation 

As explained below, the Cactus code is written in a 
modular "plug-in" way to allow for any number of for- 
mulations of the evolution system. For example, in ad- 
dition to evolving the BM system, the Cactus code has 
a straightforward ADM integrator subroutine (what in 
Cactus language we call a "thorn"), which solves the 
ADM system using a full leapfrog scheme described in 
|72| and similar to that used for evolutions in |20|] . The 
current implementation of the ADM system assumes a 
zero shift vector, and can perform conformal differenc- 
ing, as described below. We use this independent code for 
comparisons between the BM system and the ADM sys- 
tem. In this way all code infrastructure used to generate 
results is the same; only the formulation of the equations 
differs, permitting a clean comparison of results. 

The standard ADM equations are p8[: 



-2aKij + + Djfii 



(18) 



dtKi. 



—DiDja 

+a [ - {i) R l3 + Ri : 
+(3 k D k K l3 + K lk Dfi k + K kj D t (3 k 



trKK lJ -2K lk K k J 



(19) 



Here Rij is the Ricci tensor, R the scalar curvature, 
and Di the covariant derivative associated with three- 
dimensional metric g^. Note that these equations look 
much simpler than the BM Eqs. presented above, but this 
is deceptive, as the expansion of the Ricci tensor and the 
covariant derivatives brings a large number of terms al- 
ready expanded in the BM system. In fact, apart from 
the fact that the BM system introduces the V k to achieve 
hyperbolicity, the BM and ADM systems only differ by 
the introduction of first order quantities and by the use of 
flux conservative form. It is useful to notice that substitu- 
tion of the definition of V k (Eq. (Q)) into all the fluxes and 
sources detailed above allows a flux-conservative, but not 
necessarily diagonalizable, treatment of the ADM system 
as a first order system. 
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C. The Constraints 

The 3+1 decomposition of the Einstein equations re- 
sult in the evolution equations, Eqs. ( |l8|) and (|l9|), and 
additional constraint equations. These are the energy or 
hamiltonian constraint, 

R + (ti-K) 2 - K ij Kij = 2a 2 G 00 , (20) 

and the momentum constraint, 

D k (K i '-gVi X K)=aGP k , (21) 

both written here in this standard ADM form. 

Using the BM variables, we can write a more natu- 
ral way to measure the constraints for the BM formula- 
tion. The Ricci scalar term in the hamiltonian constraint 
(Eq.(pOj)) can be computed using 

r = -2d k v k + nrr,., - D k r r D k l (22) 

The treatment of momentum constraint is more subtle. 
In generating the equation for the evolution of Vi in the 
BM formulation, the momentum constraint Eq. (||l]) is 
factored in. Thus, the algebraic constraint Eq. (|J) mea- 
sures the time integral of momentum constraint viola- 
tion, since dtVi — dt(D r r i — Di r r ) = the momentum con- 
straint. Therefore, rather than measure the momentum 
constraint directly, we measure the algebraic constraint 
Eq. (||) in its place. 

D. Coordinate Systems 

We choose a 3D cartesian coordinate system with a 
general metric, general extrinsic curvature, and an arbi- 
trary 3D shift vector. In this way any slicing or shift con- 
dition may be imposed as needed. The use of cartesian 
coordinates avoids the introduction of any coordinate sin- 
gularities, and enables the treatment of many problems 
in 3D, regardless of their geometry. 

We also allow for a (time independent) conformal 
rescaling of the three-metric, which can be useful in in- 
creasing accuracy in spacetimes where the conformal fac- 
tor is known analytically, or perhaps numerically through 
a solution of the constraint equations |l9|| . The key point 
is that the derivatives of the conformal factor, provided 
in the initial data, can be known with much greater ac- 
curacy than is achieved via finite differencing on the grid 
used for evolution, and exploiting this knowledge can im- 
prove the accuracy of the evolution. In this case we write 
the metric as 

9ij=ip 4 '9ij- (23) 

This leads to a relationship between the physical vari- 
ables, denoted only here with a hat (i.e. and con- 
formal variables, 



Dki3=4\Dkij+2 J jr-gii) (24) 

Vi = Vi+4-tf. 25 

W 

We use these relationships to move the conformal fac- 
tor and its derivatives out of the flux terms and into 
source terms, allowing us to evolve the system without 
having to take numerical derivatives of the conformal fac- 
tor, while still maintaining a first order flux conservative 
form. The complete transformed equations are given in 
Appendix A. The usage of conformal rescaling is an op- 
tional parameter in Cactus, and we only use it in the 
black hole spacetime tests of Sec. |v| . 

E. Gauge Choices 

Buried in the above system of equations is the slicing 
condition. Normally considered as a supplemental condi- 
tion in the ADM evolution system, it is an integral part 
of the evolution system here, which for clarity we repeat 
here: 

d t a = -a 2 Q + a[3 r A r . (26) 

It is important to realize that one does not need to use 
this evolution equation for the lapse, as the BM formu- 
lation as presented above allows any arbitrary choice of 
lapse and shift. In principle, if one is not concerned about 
the hyperbolicity of the system, it is possible to use any 
choice and even dynamical choices that involve depen- 
dencies on the spacetime metric or the extrinsic curva- 
ture are allowed. However, given that in the future we 
are particularly interested in exploiting the hyperbolicity 
of the system, we will concentrate our studies in the the 
family of slicings introduced in Ref. ||[||. Namely, we 
admit lapses with the following gauge source function: 

Q = f(a)trK, (27) 

where the most common choices for / will be the follow- 
ing: / = 0, which implies geodesic slicing, / = 1, which 
implies harmonic slicing, f = 1/a, which gives rise to the 
so-called "1 -flog(g)" slicing. As discussed in Ref. Q, all 
choices with / > are singularity avoiding and permit a 
hyperbolic system. 

Recent work has shown the potential danger 

of hyperbolic gauges in numerical relativity, as blow-up 
along characteristics may occur depending on the choices 
for the initial data and gauge condition. This occurs in- 
dependently of the formulation of the equations. It is 
even possible in simple electrodynamics with a nonlinear 
choice of gauge. More research is necessary to character- 
ize the initial data and gauge choices that are "safe" from 
gauge pathologies. Until then, the time-honored usage of 
elliptic conditions remains the safest alternative. Maxi- 
mal slicing corresponds to the limit of divergent /. We 
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implement it in our code by not evolving Eq. Q2q), but 
rather by setting / to zero through the update step, and 
solving the elliptic gauge condition 



Aa 



K lJ K l3 a 



(28) 



after the update stage. The variables Ak related to the 
derivatives of the lapse are then computed using centered 
finite differencing. 

We also allow a non-zero shift vector. The choice of 
appropriate shift vector in 3D is still an open research 
area, and so here we demonstrate simple tests of the shift 
terms, but we do not use the shift to enforce any physi- 
cally motivated coordinate conditions (e.g., minimal dis- 
tortion). We will treat the shift as a "given" arbitrary 
function of spacetime whose derivatives are known at all 
time, which we instantaneously update every At. See 
Ref. for a full discussion of general shifts and special 
subtleties in their implementation for hyperbolic formu- 
lations (See also [ p4|j73|l for discussions of treatment of a 
general shift in another hyperbolic framework.) 



III. NUMERICAL AND COMPUTATIONAL 
CONCEPTS 

A. The Cactus Code and Computational Science 

As well as solving the Einstein equations, the Cactus 
code endeavors to address several difficult problems in 
computational science. Although these are addressed 
in detail elsewhere |pJ, [7l|,K0| , we review the basic ideas 
briefly here. 

Cactus is a parallel code, and is parallelized using the 
standard MPI message passing interface [ fr4| . This al- 
lows high performance portable parallelism using a dis- 
tributed memory model. All major high performance 
parallel architectures, including the SGI/Cray Origin 
2000, SGI/Cray T3E, HP/Convex Exemplar and IBM 
SP-2 support this programming model. Moreover, us- 
ing MPI allows computing on clusters of workstations 
using any of several free implementations of MPI. The 
parallelism software we developed in Cactus, described 
in [fn|, is a generic domain decomposition package for 
distributing uniform grid functions on various processors 
and providing ghost-zone based communications with a 
variety of stencil widths and grid staggerings. The code 
can also compile without MPI, allowing one source code 
to be used for single processor workstation development 
and for massively parallel high performance computing 
simulations. Our parallelism software is similar in spirit 
to Parashar and Browne's DAGH system |75 76 1, with 
the crucial difference being that it does not support fully 
adaptive meshes, and therefore has a much lower de- 
gree of computational complexity. However, the system 
does support the creation of multiple grids which are dis- 
tributed across all processors. This feature is used to 
provide automatic convergence testing, the importance 



of which is stressed below. The support of multiple grid 
hierarchies also allows multigrid solvers and fixed mesh 
refinement hyperbolic solvers to be built upon this paral- 
lel software. We are presently collaborating with several 
groups and colleagues to implement this and many other 
computational features which will be reported elsewhere. 

The implementation of the Bona Masso and ADM evo- 
lution equations in Cactus has been strongly optimized 
for high single processor performance on cache based ar- 
chitectures. The code very effectively utilizes the many- 
tiered memory structures of modern high performance 
computing architectures, through a variety of techniques 
described in Ref. j7fj . The combination of portable par- 
allelism with high single processor performance has led 
to a very well performing code. In recent performance 
studies, the Cactus code evolution system attained better 
than 66 GFlop/s performance on a 512 processor T3E- 
900, experiencing a speedup of more than 500 fold over 
1 processor on the 512 processor system. 

In addition to this performance related technology, the 
Cactus code attempts to be a usable code in a collabora- 
tive setting. The code has a clearly defined "plug-in" cod- 
ing style, by which users developing code to extend cac- 
tus do not modify the central code, but rather place their 
subroutines in a "thorn" which has a well defined calling 
structure. There are several positive benefits to this soft- 
ware engineering decision, as managing and maintaining 
the code becomes a distributed task. Each "thorn" and 
the central code are managed as separate modules using 
versioning software, and each small chunk has a clearly 
defined maintainer. Experimentation by a user will not 
disrupt the work of all other users, since other users will 
not be required to use new and unstable "thorns." With 
the thorn system, we are able to maintain a single cen- 
tral version of the Cactus code which all users of the code 
extend in a non-intrusive manner. 



B. Boundary Conditions 

As discussed in the introduction, boundary conditions 
are a major open research problem in numerical relativ- 
ity. It is beyond the scope of this paper to formulate an 
adequately general outgoing boundary condition. We opt 
here to use very simple boundary conditions and concen- 
trate on our evolution in the interior. We will demon- 
strate that, although poor boundary conditions can lead 
to loss of convergence in the interior of any numerically 
generated spacetime, one can still find accurate solutions 
to the Einstein equations for a finite time. In the worst 
scenario, the interior solution should always be valid for 
the Cauchy domain of dependence shown in Fig. 0, but 
generally one fares better than this with reasonable con- 
ditions, such as those we use. 
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FIG. 1. Even with an inaccurate, but stable, outer bound- 
ary, the region which is causally disconnected from the bound- 
ary has finite size after some evolution. We note that in gen- 
eral relativity the coordinate speed of propagation is not one, 
so it is possible that causal connection to the boundaries does 
not move at a linear speed, as indicated by the curved edges 
of the boundary region here. 

The boundary condition we use is a simple copying 
boundary (sometimes called zero order extrapolation). 
That is, for each point on the physical outer boundary, 
we copy all the variables from the point nearest the in- 
side. In practice, this condition will prove very effective 
in several scenarios. It has the effect of canceling the flux 
difference in the exterior (as the finite differences of the 
last points will be zero). This is a valid approximation 
to "outgoing" boundary conditions when the boundary 
is close to linear perturbation around flat spacetime. In 
this special case, the sources of the BM system are close 
to zero and the system approximates a set of linear wave 
equations, so canceling the exterior flux effectively pre- 
vents any incoming information from outside the domain. 

Following Ref. pH , we allow octant boundary condi- 
tions, which are appropriate for spacetimes with rota- 
tional and equatorial plane symmetry. This allows us to 
simulate black hole spacetimes with symmetries as full 
3D problems, while using one eighth the computational 
resources necessary when evolving on a full grid. Many 
interesting problems, including Schwarzschild [ pl| , ax- 
isymmetric black hole collisions [T7| and distorted black 
holes j2J], and even some full 3D data sets with cer- 
tain dependence on the azimuthal angle [ p3|j78| , can be 
treated with this symmetry, allowing a great savings in 
computational resources. Of course, our code can run 
without this boundary condition also, and as demon- 
strated through comparisons running full and octant 
grids the use of octant symmetry does not affect 
results. 

There are many other boundary conditions which 
are applicable to three-dimensional numerical relativity, 
which we do not consider here. However, they are worth 
mentioning. The apparent horizon boundary condition 
(AHBC) (3^^] adds a boundary at the causal interior 
of a black hole spacetime. Recent progress on the outer 
boundary treatment, such as matching schemes to per- 
turbative B or characteristic 1791 evolution schemes, look 



very promising, and could be ultimately used by Cactus. 
Another promising boundary treatment involves mov- 
ing the outer boundary to infinity pfj| by conformally 
rescaling the metric, as per Friedrich's hyperbolic system 
J49| . This has proven very successful in one-dimensional 
calculations |H],^2| and higher dimensional calculations 
with this method should be available soon. Finally, we 
reiterate that boundary conditions are a major motiva- 
tion for hyperbolic treatments of the Einstein equations. 
Through study of the eigenfields and eigenvalues of the 
transport system they provide more information about 
the flow of information at the boundaries, which can be 
exploited in numerical methods pi] ]. 

The interior of black holes is usually handled with an 
isometry condition, which identifies the interior of a black 
hole with the isometric exterior via inversion through the 
sphere. This has been crucial in numerous black hole 
evolutions published to date (see, e.g., @|20g[r§]. 
We do not use a three-dimensional isometry condition, 
as is described in |^l|,^0|, since we must transform not 
only the metric and curvature tensor, but also the first 
order quantities (such as derivatives of the metric and 
the vector Vk) which do not transform as a tensor. An 
isometry could be implemented for the BM formulation 
in principle, but as we are looking to move to more gen- 
eral methods to treat the black hole, such as an AHBC, 
we have chosen not to do so at present. Furthermore, 
as shown in Ref. [jilt , with certain slicing conditions like 
maximal slicing, even without AHBC the isometry con- 
dition can be ignored and both regions inside and outside 
the horizon can be evolved, as long as the lapse collapses 
sufficiently quickly in the vicinity of the singularity. We 
will make use of this property of maximal slicing in tests 
presented below. Recent work proposes an alternative to 
isometry conditions by "stuffing" with matter the inte- 
rior of black holes ( "stuffed black holes" ) and we are 
currently investigating this approach for 3D spacetimes. 



C. Evolution Schemes 

1. The Strang Splitting 

Following the numerical discussion of the BM system in 
Ref. we will split Eq. (|l|) into two separate processes. 
The transport part is given by the flux terms 



d k F k u = 



(29) 



The source contribution is given by the following system 
of ordinary differential equations 



d t u = S-U 



(30) 



Numerically, this splitting is performed by a combination 
of both flux and source operators. Denoting by E(At) 
the numerical evolution operator for system ([!]) in a sin- 
gle timestep, we implement the following combination 
sequence of subevolution steps: 
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E{At) = S(At/2) T(At) S(At/2) (31) 



3. The Flux Update Methods 



where T, S are the numerical evolution operators for 
systems (|29| ) and (]3Cj), respectively. This is known as 
"Strang splitting" |72| . As long as both operators T and 
S are second order accurate in At, the overall step of 
operator E is also second order accurate in time. 

This choice of splitting allows easy implementation of 
different numerical treatments of the principal part of 
the system without having to worry about the sources 
of the equations. Additionally, there are numerous com- 
putational advantages to this technique, as discussed in 
frO| . Theoretical and practical advantages for general 
relativistic hydrodynamics, where the source step cou- 
ples the equations for the whole system of Einstein plus 
matter equations, will be detailed elsewhere. 



2. The Source update method 

Currently, we treat the source integration with a sec- 
ond order predictor-corrector method ]7^] . During this 
step, we only need to evolve the 16 quantities which have 
a source (gij,Kij,V k and a). 

We use standard finite difference notation here. Sub- 
scripts denote grid index, and superscripts denote time 
index. For instance, u™ k is the value of field u at spa- 
tial grid point i, j, k and time level n. We use the special 
upper indices p and c to denote the predicted and cor- 
rected values during an update cycle, as we define below. 
In order to update a variable u (running through the 16 
quantities with source) at time level n to the future time 
level n + 1, we first compute the "predicted value" u P j k 



at every point i, j, k of our computational grid. 



(32) 



where 



is u at current time step n and grid point 



i,j,k, and At is the time discretization interval. With 
this predicted value of u p , we compute the predicted 
sources and take a corrector step: 



\j,k 



AtS(ul 



3,k 



)• 



(33) 



Finally, the evolved value of u at the next time step n + 1 
is the average of the value at time step n and the correc- 
tion: 



,n+l 

%j,k 



i,j,k 



)/2 



(34) 



In practice, the steps (|33| ) and ([34]) can be combined into 
one. Note that this is a completely local operation at 
every grid point, which allows a high degree of optimiza- 
tion 170]. Higher order methods for source integration 
can be easily implemented, but this will not improve the 
overall order of accuracy. However, in special cases where 
the evolution is largely source driven p6[ , it may be im- 
portant to use higher order source operators, and this 
method allows such generalizations. 



The implementation of numerical methods for the flux 
operator is much more involved, and we have many 
choices at our disposal, ranging from standard choices 
to advanced shock capturing methods [ pl| , ^5p |. In this 
paper, we will limit ourselves to two methods: the Mac- 
Cormack method, which has proven to be very robust in 
the computational fluid dynamics field (see, e.g., Ref. |||] 
and references therein), and a directionally split Lax- 
Wendroff method. These schemes are fully second order 
in space and time. Although the Cactus code has a mod- 
ular structure allowing numerous numerical methods to 
be plugged in and applied to problems for which they 
may be best suited, in this paper we restrict ourselves to 
results with these two methods. Unless otherwise noted, 
results are generated with the MacCormack method; use 
of the Lax-Wendroff solver will be explicitly noted. 

Following the previous notation we define our fluxes in 
individual directions x, y, and z as FX, FY, and FZ 
respectively. 

The MacCormack method evolves a given quantity u, 
which now runs through the 30 dynamical variables hav- 
ing fluxes {Kij,D ki] ,V k ,A k ; the A k and D kij only have 
fluxes in one direction, which is explicitly exploited in 
our code) with the following algorithm: First, in order to 
update the variable u to the time level n + 1, we compute 



with first order backward fi- 



the "predicted value" u; 
nite differences: 

u>ij,k = <j,k + ^( FX «j,k) ~ FX W-i,j,k)) 

+ ^ FY « 1 ,k)- FY « ] -i,k)) (35) 

+ ^(FZ(ul^)-FZ(u^ k _ 1 )) 

where, in addition to the quantities defined above, Aa;, 
Ay, and Az are the spatial discretization intervals. Note 
that this predicted step can be done in a given direction 
(say x), from grid points 2 to nx (total number of grid 
points in that direction), as the first order backward dif- 
ferencing only requires i — With this predicted value of 
u p , we recompute predicted fluxes and sources and take 
a corrector step with forward finite differencing: 



At 



l i,j,k ~ u iJ 



k + Xz( FX (< + i, 3 ,k)-FX(u^ k )) 



Ax 

+ ^(FY(ul j+ljk )-FY(ul jik )) (36) 
+ ^ { FZ(ul j<k+1 )-FZ(ul jtk )) 

Now we can correct the interior points of the domain 
from 2 to nx — 1, as we have a prediction for the last 
plane at nx. Finally, the evolved value of u at the next 
time step n + 1 is the average of the value at time step n 
and the correction: 
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,n+l 



(«?j,fc+«Ufc)/ 2 - 



(37) 



A similar method could be obtained interchanging the 
order or backward and forward derivatives in the pre- 
dictor and corrector steps. We note that both methods 
can introduce certain spatial asymmetries in a numerical 
evolution, due to the preferred order of finite difference 
operations in the predictor and corrector steps. These 
asymmetries converge away to second order, as we will 
discuss below when presenting results. 

The directionally split Lax-Wendroff method uses a 
series of one dimensional Lax-Wendroff integrations to 
complete a full three dimensional integration step. In 
one dimension, the Lax-Wendroff scheme is 



v 



n+l/2 
-1/2 



1 



u? +1 = U? 



Ax y y 1 



At 
2Ax 

) 



(F(u? +1 ) - F(u?)) (38) 



fl/2 
-1/2 



F(u' 



1+1/2 
-1/2 



))• 



(39) 



Several options exist to turn Lax-Wendroff into a 
three dimensional scheme. Here we choose direc- 
tional splitting ||72t . Defining X(At) to be a one di- 
mensional Lax-Wendroff in the x— direction, Y(At) in 
the y— direction and Z(At) in the z— direction we de- 
fine a full flux time step as X (At)Y (At) Z (At) on 
the first step, Y(At)Z(At)X(At) on the second step, 
Z(At)X(At)Y(At) on the third step, and then repeat 
the prescription. This permutation leads empirically to 
a second order in space and time scheme, as we shall 
demonstrate below. The advantage of this directionally 
split Lax-Wendroff is that, by turning the problem into 
a set of one dimensional PDEs, implementation of a sim- 
ple inner (apparent horizon) boundary condition becomes 
easier, as will be reported elsewhere |87 . 



D. Convergence 

Since the pioneering work of Choptuik [Q, the usage 
of convergence tests in numerical relativity is slowly be- 
coming standard practice fl41^|^,f45|] . The recent discov- 
ery and characterization of gauge pathologies j|5| stresses 
the importance of careful convergence analysis, especially 
in 3D numerical relativity, as simulations may hide solu- 
tions that "look" reasonable for a given resolution but 
do not satisfy the constraints. For completeness, here we 
review the basis of convergence tests. We will discuss the 
case of numerical discretization of PDE's with finite dif- 
ferences. Similar arguments can be developed for other 
approaches. 

Assuming that we have well behaved solutions which 
allow a expansion in Taylor series, we can relate the nu- 
merical solution S to the analytical solution S in the 
following way: 



S = S + 0(A° 



(40) 



where A is the grid spacing. Consistent numerical simu- 
lations must demonstrate that some form of this relation 



is obeyed, as the refinement of the grid should always 
improve the solution. In many cases, it is actually pos- 
sible to measure the convergence rate a. This analysis 
is crucial if one is to understand how close a given nu- 
merical solution is to the true analytic solution, which is 
generally not known. 

Given three discretized solutions, S'(A), S(A/q) and 
S(A/q 2 ) we find that 

L = S(A/q)-S(A) = 0((A/qr - A°) (41) 
M ee S(A/q 2 ) - S(A/q) = O ((A/q 2 Y - (A/qf) . (42) 

We define precisely the intuitively clear "-" operator be- 
low. Dividing and canceling A° we find 



L q- a - 1 



M q- 2a - q-° 
so solving for a, 

logg 



(43) 



(44) 



Eq. (|4|) is the principal definition of the convergence 
rate a that we will use. In practice, we use q = 2 for our 
convergence tests, that is, we double or halve our grid 
resolution in a sequence of simulations when determining 
a. 

The definition of the "-" operator used to form L and 
M is very important. If the points of S(A/q) and S(A) 
are coincident on the A grid, then simply pointwise sub- 
traction followed by a norm of the difference can generate 
the "-" operator. We can schematically represent this op- 
eration as 

S(A/q) - 5(A) = \S(A/q) qhqMk - S(A) itj , k \, (45) 

although often more complicated index juggling that sim- 
ply i — > qi is required. We note that \x\ denotes some 
norm over the k space (e.g., maximum, C\, £2)- 

If the points are not coincident, a possibility is to use 
some norm over the solutions and then define "-" as the 
difference of those norms. That is 



S(A/q)-S(A)EE\S(A/q)\-\S(A)\. 



(46) 



We call this convergence in the norm. This method has 
the advantage that it is very easy to calculate during run- 
time of a parallel code, but is often susceptible to large 
amounts of noise. If we have an interpolation operator 
which interpolates a solution from a grid with res- 
olution A/q to one with resolution A, we can define 

S(A/q) - S(A) ee \l£ /q S(A/q) - S(A)\. (47) 

Generally, an interpolator of at least order a is required 
to do this style of convergence testing. 

Finally, when an exact solution is known, we only re- 
quire two numerical solutions to the equations to mea- 
sure a. That is, given a discretized solution at S(A) 
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and S(A/q) and an exact solution S, we can form two 
differences pointwise, 



H(Ax) = + £(A CT ), 



(51) 



L = S-S(A) = 0{A") 
M = S- S(A/q) = 0{q- a A a ) 

and therefore find the relationship 

L = q a M, 



(48) 
(49) 

(50) 



and again we recover a from Eq. (|4|). Simply said, for 
a second order method, the error should be four times 
larger on the coarser grid than the finer grid. This 
method will prove valuable for calculating convergence 
against known solutions, convergence of constraints, and 
convergence of fictitious numerical errors, such as asym- 
metries. 

As before, we are faced with the problem of comput- 
ing the quotient L/M accurately, especially in the case 
of fields which go to zero. Once again, we can solve this 
problem by interpolating M onto the grid which L in- 
habits and forming the quotient pointwise. 

When the convergence rate is expected to be second 
order, with this technique we can also measure a graphi- 
cally at all points. That is, if we have a known solution 5 1 , 
we can plot, S(A)-S, (S(2A)-5)/4, (S(4A)-S)/16 and 
so forth. If the points agree, then we have second order 
convergence. This method has the advantage that point 
to point noise present in calculating a can be eliminated 
"by eye," and we shall use this method often below. 

Convergence testing is an essential component of a 
battery of code tests. Demonstrating that an evolu- 
tion scheme has the appropriate convergence order shows 
that boundary treatments, methods, and infrastructure 
are coded properly. Studying convergence properties can 
help diagnose and track subtle errors in a code. How- 
ever, showing that, for example the metric function g xx 
converges does not imply that one is solving the Einstein 
equations; it merely means that one is solving the coded 
evolution equation to second order. Thus, convergence 
testing against known solutions is important. In a few 
rare cases, notably a geodesically sliced black hole, there 
are exact solutions to the non-linear dynamical Einstein 
equations. In this case one can show not only that nu- 
merical results converge to something (that is, we find 
a = 2 using definitions (|44|)), but also that they converge 
to the right thing (that is, we get order a = 2 when com- 
paring against the known solution, using the definition 

With the Einstein equations, however, we play on fa- 
vorable ground, as we always have an analytic solution 
at our disposal: the vanishing of the constraints. That is, 
all correct solutions to the fully nonlinear Einstein equa- 
tions have the property that the hamiltonian and mo- 
mentum constraints must identically vanish. Regardless 
of the relativistic system being simulated, if the initial 
data satisfies the constraints, then so must all subsequent 
time steps. We assume the behavior of the hamiltonian 
constraint, H, is 



where E{A a ) is the error due to finite differencing with 
a spatial step A. Choptuik has investigated this point 
at great length in Ref. |pgj| , where he shows that for a 
consistent finite differencing of the free evolution of the 
Einstein equations, the constraints have the same order 
error as the evolution scheme. Choptuik demonstrated 
this in spherical symmetry, and here we demonstrate this 
in full three dimensional numerical relativity. With rela- 
tion Eq. (fil]) in hand, forming L and M in the language 



of Eq. (49) simply amounts to looking at the value of the 
constraint. If we double the resolution, and the numeri- 
cal code is solving the Einstein equations, our constraints 
must drop by a factor of four (for a second order scheme) 
everywhere. (We note that one may use the constraints 
to eliminate one of the evolution equations, and with this 
approach it would be reasonable to expect that the code 
could demonstrate an independent construction of the 
eliminated evolution equation converging to zero, rather 
than the constraint.) 

Having discussed the convergence techniques we use to 
study the performance of the Cactus code, we describe 
briefly our philosophy of their use before moving on to 
examples below. An important point is that a 3D code 
should exhibit convergence, even if the resolution is too 
low to exhibit a high degree of accuracy. For instance, a 
given numerical result may differ from the true analytic 
solution by a large factor. This is not necessarily a major 
concern, as long as doubling the resolution can quarter 
the error. By running simulations at different resolutions, 
one can then estimate how close the numerical solution is 
to the analytic solution, understand the behavior of the 
truncation error, and estimate the resolution required to 
obtain a solution to the desired accuracy. We regard this 
as a crucial requisite of a code, which as we show below, 
our code has. 

We note that it is possible that resolution can be too 
low to allow an evolution beyond a certain point. For 
instance, a geodesically sliced black hole with very low 
resolution may crash before (or after) a higher resolution 
simulation, since there is a physical singularity present. 
However, in this case convergence should show the re- 
gion in which the simulation is accurate. When a solu- 
tion starts failing to converge, the evolution is pro bably 
about to fail. We will see an example of this in Sec. VI C . 
Convergence at the boundaries also offers useful infor- 
mation. Using our simple "copying" boundary condition 
without any advanced treatment of the system, we expect 
that condition to have a first order effect on phenomena 
which interact with the boundary. Finally, we regard sec- 
ond order convergence to be desirable, but not necessary, 
in order to verify a given numerical result. For instance, 
often one does not have second order convergence near 
boundaries. But such an effect can be studied and un- 
derstood. The key point is that one should know that 
a code converges at or above the expected order, even if 
that order is one. 
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IV. FLAT SPACE TESTS 
A. Dynamically sliced Flat space 

One of the crudest first tests of any 3D cartesian based 
code is, given geodesically sliced Minkowski space, does 
the code produce 1 and forever. Of course Cactus does, 
but this test is almost useless, since one cannot measure 
convergence, and all constraints are trivially satisfied. A 
more interesting and much more important test is that of 
a dynamically sliced flat spacetime. That is, we choose 
an initial lapse in Minkowski space which is not unity 
everywhere, and then we evolve this system with a "live" 
slicing condition on the lapse a, such as harmonic slicing. 
Such an idea has been suggested in the past by York J4q] 
and also implemented and studied in detail by Masso [p6| . 
We have examined three distinct cases, ID periodic data 
(e.g., a is a periodic function of one coordinate only), 3D 
periodic data, and 3D data where a falls off to unity at 
large radii. The first two allow us to study Cactus with- 
out boundary affects, and the last allows us to evaluate 
the quality of our boundary conditions. In this section we 
present 3D simulations with "copying" boundary condi- 
tions and harmonic slicing, and in the following sections 
we also chose various shifts in both one and three dimen- 
sions. Simulations here are performed with the Einstein 
system (n = 1). 

We note this problem is similar to the example used 
to study coordinate conditions discussed in }43|], but 
with some important differences. First, we evolve with 
harmonic slicing throughout the entire evolution, rather 
than using the lapse to generate a small "bump" in trK , 
followed by maximal slicing, as in Ref . E] . Secondly, in 
harmonic slicing of flat space, the lapse evolution equa- 
tion becomes wavelike, so our initial pulse travels off the 
grid as a wave pulse. 

In the 3D case, we choose an initial lapse with a gaus- 
sian bump specified by 



a = 1 + A exp 



(52) 



In the ID case, we choose an identical form, simply re- 
placing r with x, y, or z. In both cases, we see simple 
wave-like propagation in the lapse, as demonstrated be- 
low. For our first test, we evolve this dynamically sliced 
flat space system in 3D with a resolution Ax = Ay = 
Az = 0.01 on a grid of 101 3 centered around r = 0. 
We choose the parameters A = 0.05 and a 2 = 0.05. 
Two-dimensional slices in the z— plane of the evolution of 
this initial lapse in a harmonically sliced spacetime, with 
"copying" boundary conditions discussed in Sec. HI B 
above, are shown in Fig. ^|. Other metric functions, al- 
though initially taking a Minkowskian form, develop sim- 
ilar dynamics. 
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FIG. 2. We show slices of the evolution of the lapse in the 
3D dynamically sliced flat spacetime described in the text. 
Figures (a), (b), and (c) show slices in the z = plane at 
times t = 0.0, 0.1, and 0.2. In Figure (d) we show the value 
of a on the y = 0, z — line evolving in time as a colored 
contour map. We note that in this simulation, the majority 
of the pulse has not yet hit the boundary of the computa- 



tional domain, 
of Ax = 0.01. 



101 grid points were used with a resolution 



We demonstrate that the hamiltonian constraint con- 
verges at second order in the interior in Fig. |^, where we 
show the constraint at three different resolutions, with 
the appropriate factors of four and sixteen. The fact that 
the lines are coincident demonstrates second order con- 
vergence. The actual value of the convergence exponent 
on the grid is above 1.9 for the entire evolution, until the 
pulse interacts strongly with the boundary. 
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FIG. 3. We demonstrate convergence of the hamiltonian 
constraint to zero in the interior of the 3D evolved dynam- 
ically sliced flat space. Although the a:— axis is shown here, 
other directions have similar results. The fact that the high 
resolution hamiltonian is equal to one-quarter the medium 
resolution, and this one is itself one-quarter the lower reso- 
lution, indicates that the hamiltonian converges rigorously to 
zero at second order in the interior. However, the logarithmic 
scale reveals, at a very low level, a lack of second order con- 
vergence near the boundaries at later times. This is caused 
by our outer boundary condition, which is not expected to be 
second order accurate. 

We noted above that due to the upwind/downwind na- 
ture of the MacCormack predictor-corrector method we 
use, certain asymmetries in the evolution are introduced. 
In Fig. ^ we see that symmetry around the origin of the 
coordinate system is not maintained except in the limit of 
a converged solution. (We note rotational symmetries are 
obeyed. By this we mean that, given symmetric data, our 
code will generate identical solutions along an x-directed 
and y-directed slice of our data. However both of these 
solutions will be (identically) asymmetric around the ori- 
gin.) This asymmetry is purely an artifact of our method 
having an upwind/downwind nature, as shown in the fi- 
nite difference representation. As such, this asymmetry 
should converge away at second order. In Fig. we show 
that this asymmetry is an artifact of numerical error, and 
consequently, converges to zero by measuring the asym- 
metry, E = a(x) —a(—x), for the evolved flat space case. 
Clearly this should be zero in the converged limit, so the 
numerical solution should obey E(Ax/2) — E(/S.x)jA if 
our method converges at second order. From Fig. we 
see that this relationship is obeyed except at the bound- 
aries, where our boundary condition imposes a first order 
asymmetry on the system at late times. 
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FIG. 4. We show symmetry violation due to the MacCor- 
mack predictor corrector method, and how that converges 
away. For the dynamically sliced flat space model, we show 
E = a(x) — a(—x) for x > at the times t = 0.8, t — 0.16 
and t = 0.24. We show E/i at a low resolution (dotted line) 
and E at twice the resolution (solid line). The fact that the 
high resolution error is less than or equal one-quarter the low 
resolution error indicates that the method's asymmetry con- 
verges to zero at second order. An interesting feature of this 
figure is that it demonstrates the first order nature of our 
boundary condition clearly. Since the lapse becomes dynamic 
on the boundary at later times, convergence order drops from 
two (which the evolution system obeys) to one (which the 
boundary condition obeys) as the wave propagates towards 
and through the boundary. 

Figs. H and [| also give an interesting indication of our 
boundary conditions when dynamics are present at the 
boundaries. As shown in Fig. ||, the traveling pulse in 
the lapse is approaching the boundary by late times in 
our simulation. Once the dynamics reach the boundary, 
convergence drops from second order towards first order 
there. This is indicated in Fig. [| by the high resolution 
case (solid) having more than one-quarter the error of 
the low resolution case (dotted line), and by non-second- 
order convergent (although small) errors in the hamilto- 
nian constraint, as shown in Fig. |3|. That is, the solid 
line is above the dotted line. 

So far, we have only measured convergence of metric 
functions and constraints. We can also examine the phys- 
ical properties of our underlying spacctime. In this space- 
time, we can demonstrate that we are evolving Minkowski 
space by measuring the Riemann invariants I and J, com- 
puted using a 3 + 1 method |8!| . These should be identi- 
cally zero, but they will not be due to finite differencing 
errors. However, we can test how they behave with vary- 
ing resolution. In Fig. || we show |/| at three different 
resolutions for the distorted flat space case considered 
here. We note that, firstly, \I\ is small, and also that it 
decreases faster than second order with grid resolution 
towards zero. In fact, in this case the convergence expo- 
nent for |/| is very close to four. Clearly boundary effects 
are evident, driving the system away from the underlying 
flat space. 
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FIG. 5. We show the Riemann invariant \I\ along the a;- line 
in the dynamical flat spacetime at various resolutions. We 
show the invariant at three resolutions with appropriate fac- 
tors of 1/4 and 1/16, and note that the invariant is converging 
towards zero faster than second order. In this case, the con- 
vergence exponent is very close to four. 



FIG. 6. We show propagation of the lapse in a dynamically 
(harmonically) sliced Minkowski space where the initial lapse 
is chosen to be a periodic gaussian function of z alone. We 
see the expected affect of the shift. That is, shifts of ±1 force 
the propagation of the lapse to have zero coordinate speed in 
one direction. 



B. Testing the Shift 

We now introduce the shift vector j3 l to test its effect on 
the solution. The dynamically sliced flat spacetime is an 
excellent case to test the shift terms in Cactus. We first 
examine a constant shift, and then move to a spatially 
varying shift to test all terms related to the shift vector. 

1 . A Test of a Constant Shift 

As a first simple test, we chose the one dimensional 
periodic initial lapse, and evolve this on an explicitly ID 
grid with periodic boundary conditions (that is, we use 
Cactus on a (nx, 1,1), (l,ny, 1) or (1,1, nz) sized grid). 
The initial lapse is chosen the same as in Eq. (^2|), with 
r replaced with x, y, or z alone, with a 1 = 0.05. In 
this harmonically sliced system with a constant shift, the 
evolution equations become wavelike for the lapse, with 
the propagation velocity being 1 ± (3. 

In Fig. g we see exactly this propagative behavior. The 
lapse function a is shown for three cases, (3 = and /3± 1. 
For f3 — 0, the wave propagates with speed c = 1 in both 
directions. For the shift chosen as ±1 we see the speed of 
the waves to be two or zero, depending on the direction 
of propagation and the sign of the shift. This can be 
clearly read from the graph, where the propagation in 
the t direction (vertically) is 0.5 in all cases, and the 
propagation distance in the z direction is 0.5 in the zero 
shift case, and 1.0 and 0.0 in the ±1 shift case. The other 
metric functions, not shown, exhibit similar behavior. 



2. A Test of a Spatially Dependent Shift and an Important 
Lesson 

We next turn to a spatially non constant shift as a test 
of our code, 

[3 X = py = [3 z = Ae- {x2+v * +z2)/a2 . (53) 

We here only consider the cases of A < 1, a sub-tachyonic 
shift. The gaussian width s is chosen so the shift is 
resolved but effectively vanishes before the boundaries. 
This choice of shift will test all terms in our (non- 
conformal) evolution equations, since it has derivatives 
of all shift terms in all directions. The following runs 
were performed with Ax = Ay = Az = 0.01, a 2 = 0.02, 
A = 0.5, and with 101 grid zones in each direction. 

Using this shift, we discovered an error in our code, 
which is worth discussing. In an initial version of our 
code, we had an error in the shift term for the sources of 
the V variables. Rather than the correct term, 

2{D ri s - 6?D^ r )B r a (54) 

we had the different, although very similar, 

■1:1),;- <<;/>•'..,) />'.'• (55) 

(Recall that B is not symmetric.) As we show now, by 
only performing convergence tests we were able to diag- 
nose and track down the code error, without appealing to 
any analytic solutions beyond the vanishing of the con- 
straints. 
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In Fig. [?] we show the evolution after some time choos- 
ing the shift in Eq. (|53|), with and without the error 
above. As is clear, the evolutions are very similar; in 
fact, had two different codes given this result, without 
further testing one would be tempted to say the results 
are the "same" and so the codes "agree" . 
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FIG. 7. We show two evolutions of a in our distorted flat 
space model with a spatially dependent shift, using the Ein- 
stein equations in one case, and the equations with a small 
error in the second. In (a) we show the numerical solution 
after eight iterations for the case with the correct shift terms 
with a solid line, and the results with an error in the shift with 
a dashed line. At this level, the plots are indistinguishable. 
In (b) we show the difference between the two evolutions, and 
notice the difference is negligible compared to the disturbance 
in the lapse. 
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FIG. 8. We show the convergence of the hamilto- 
nian constraint for the Einstein equations above and the 
non-convergence of the constraint for the Einstein equations 
with an error below. We note that even though the error 
in our lapse evolution is very small, the convergence simply 
fails for the incorrect equation (note in the lower plot that the 
hamiltonian is the same for both resolutions, although the fig- 
ure might mislead the reader because we introduce the factor 
of 4 that we would expect for convergent results). Again, this 
demonstrates second order convergence for the correct equa- 
tions. We note that in both cases the hamiltonian constraint 
is "large"; about 0.2 in the high resolution correct case (0.8 in 
the low resolution case) and about 15 in the incorrect case (be- 
ing non convergent, stays the same for both grid resolutions). 
The only way to determine if the constraint is too "large" is 
to test its convergence towards zero, which is a feature of only 
the Einstein equations in this case. 



However, in Fig. |8| we show that the hamiltonian con- 
straint, as defined by Eq. (^|]), converges to zero for the 
Einstein equations, and fails to do so for the system which 
is not. The failure to converge is clear and large. We 
note that even with fairly low resolution we can demon- 
strate that our code is correct or incorrect by showing 
merely the convergence of the constraints and we did not 
need an exact solution for the spacetime (other than the 
vanishing of the constraints). We feel that this clearly 
demonstrates that convergence testing constraints is an 
important and strong test of any code. 



V. WAVE SPACETIME TESTS 

Although hyperbolic reformulations of the 3D Einstein 
equations have not been used in a wide variety of space- 
times before this publication, they have been applied to 
linearized wave spacetimes |4l],|7j . The current version of 
this code owes much to the implementation of the "H" 
code described in Ref. |x]]. As we reviewed in the in- 
troduction, this "H" code used a previous BM formula- 
tion of the equations that required the exclusive use of 
harmonic slicing and zero shift vector O. That code is 
now obsolete, although all the tests of the "H" code de- 
scribed in Ref. can be replicated successfully by this 
new and much more advanced version of the code. All 
the tests presented here are run with the "Ricci" system 
(n = 0), as this corresponds more closely to the sim- 
ulations performed with the "H" code. Here we detail 
some of these comparisons, evolving linear initial data 
that describe weak gravitational waves. The interesting 
transition from linear to non-linear effects described in 
Ref. Q will not be studied here, although it is possible 
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to reproduce those effects with the two formulations (BM 
and ADM) implemented in Cactus. 

Further studies of stronger gravitational wave inter- 
actions and their possible collapse to a black hole are 
underway and will be described in a future publication 
in this series, where appropriate slicing conditions for 
wave spacetimes will be considered in detail. In this sec- 
tion, we will focus on two cases, colliding plane waves 
and quadrupolar waves, and limit our gauge to harmonic 
slicing. 

A. Plane Waves 

We consider linearized plane wave solutions, following 
the test in section III of Ref. fji] . The line element is 
written 

ds 2 = -dt 2 + (1 + f(t, z))dx 2 + (1 - f(t, z))dy 2 + dz 2 . 

(56) 

For small /, the linearized Hamiltonian constraint is sat- 
isfied, and the evolution of the spacetime is governed by 
the linear wave equation 

d?f(t,z)=d 2 J{t,z), (57) 

that describes plane waves propagating in the z direction. 
We use the Gaussian-shaped packet: 

f(t, z) = A Re -(Mt-*-»)M* cos f^ {z ^ 

+ A L e-^ t+z - a ^^ 2 cos (^-(z + t)^ , (58) 

The amplitudes An and Al represent the amplitudes of 
waves traveling to the right and left, respectively, with 
a Gaussian shape of width a and centered at z = ±a at 
t = 0. A is the wavelength of the Gaussian-modulated 
oscillations. 

In Fig. |^ we show the evolution of the metric compo- 
nent g xx for a single wave moving in the —z direction. 
We have chosen the shape parameters a — 2, A = 1,Al — 
0.00001, A R = and a = 3, with Az = 0.025. This figure 
replicates Fig. 1(a) of Ref. which used an ADM code 
with an staggered-leapfrog algorithm. We notice that the 
wave is transported with a small loss of amplitude, due 
to dispersive effects in the MacCormack predictor cor- 
rector scheme. The measured convergence a is very close 
to 2. This is one of the simplest tests of a numerical 
code designed to evolve waves and the results obtained 
are in agreement with the well-known numerical proper- 
ties of our standard methods applied to the linear wave 
equation. 
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FIG. 9. We show the evolution of g xx for a single plane 
wave moving in the — z direction at times t = 0, t = 3 and 
t = 6 for a gaussian wave packet. This figure replicates 
Fig. 1(a) of Ref. ^J. The dispersive nature of the Mac- 
Cormack method can be appreciated in the non-symmetric 
propagation of the gaussian packet. 

A more involved test results from colliding plane waves. 
Unlike the previous test, in this case we deal with non- 
trivial spacetimes: theoretically, it is known that such 
spacetimes will develo p a singularity in the future (in 
the non- linear regime) |9l]j92]| ; numerically, coupled non- 
linear and finite differencing effects can lead to spurious 
numerical evolution |4l| . Hence, they provide a stronger 
test of a numerical code. In Fig. [l0| we show the evolu- 
tion of a colliding wave system. Two wave packets orig- 
inally start, moving inwards, centered at z = ±3. We 
choose the same parameters as the single wave packet 
except for the amplitudes Ar — Al — 0.025. The pack- 
ets collide at the center at time t — 3 and then continue 
on. Once again, dispersion is visible when the waves 
return to their original images at t — 6. This figure repli- 
cates Fig. 6(d) of Ref. (4^]. There it was shown that the 
staggered-leapfrog method was prone to a large secular 
drifting after the packets collided, which does not occur 
with our MacCormack method. 
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FIG. 10. We show the evolution of g xx for two colliding 
plane wave packets. At t — the two packets are centered at 
z = ±3, they collide and superimpose at t = 3. At t = 6 the 
left and right packets have interchanged their positions and 
should be coincident with their shapes at the initial time. The 
difference is due to numerical dispersion. Continuing the evo- 
lution, the packets are more dispersed at t = 9. This figure 
replicates Fig. 6(d) of Ref. jy]. The MacCormack method 
used here does not exhibit the drifting after the collision ex- 
hibited by the staggered-leapfrog method in Fig. 6(a) of that 
reference. 



B. Pure Quadrupolar waves 

The numerical simulation of quadrupolar linearized 
wave solutions to the Einstein equations has been estab- 
lished as an standard test of 3D numerical codes |?],|d],[)3| . 
One of the reasons of their appeal is the existence of a 
family of analytic solutions for both even- and odd-parity 
and the independent azimuthal modes fl94| . But more 
importantly, we also need to model their evolution ac- 
curately, as quadrupolar modes are a dominant signal in 
the late time evolution of black hole spacetimes. In this 
section we compare evolutions of quadrupolar waves in 
Cactus with previous results, following again the exten- 
sive tests and discussions of Ref. [[ll| . Due to the length 
of the analytical expressions, we do not write the solu- 
tions here and refer to the reader to Ref. (94| or Ref. 

©■ 

We start by evolving even-parity waves with an am- 
plitude of 10~ 5 and quadrupole numbers I — 2 and 
m = 0. The details of this setup are given in section 
VI of Ref. Q. In Fig. |l] we show the evolution of g X x 
along the z-axis performed on a grid of 120 3 points with 
Aa; = Ay = Az = 0.05. This replicates Fig. 9(c) of Ref. 
EjJ . We can see how an initially moderate wave packet 
near the center of the grid oscillates and propagates off 
the grid, as expected. 



FIG. 11. The evolution of metric function g xx along the z 
line is shown for linear quadrupolar waves with I = 2, m = 
and a low amplitude packet, which corresponds to a pertur- 
bation of 0.025% in the metric functions. The wave expands 
outward as time increases, returning to a flat profile after 
t = 4. This replicates Fig.9(c) of Ref. @. 
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FIG. 12. We show the time evolution of the wave-like quan- 
tity r{g xx — 1) measured at the outer boundary for the sim- 
ulation shown in the previous figure. The wave pulse arrives 
to the boundary at around t = 4, oscillates and leaves the 
computational grid. This serves as indicator of the outgoing 
condition provided by our simple copying boundaries 

In Fig. [l^ we show the time evolution of the quantity 
r{g X x — 1) at the outer boundary of our grid, as an in- 
dicator of the clean outgoing condition provided by our 
simple copying boundaries. This measure of the wave 
simply separates the perturbation from the background 
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Minkowski metric, and corrects for the 1/r falloff. It 
is not a gauge-invariant measure of waves, such as that 
used in Ref. 78 1. Detailed studies extending these re- 
sults beyond linear wave regimes are under way and will 
be published elsewhere. 

Ironically, the extensive work of Ref. jyj does not in- 
clude results with any truly 3D spacetime, as the cases 
studied for quadrupolar waves correspond to axisymmet- 
ric waves of azimuthal number m = 0. In this paper 
we will extend the results of that reference by setting 
up a slightly more realistic scenario, tuning the parame- 
ters to mimic what we expect from late time ringdown of 
black hole simulations. Therefore, we will evolve non- 
axisymmctric quadrupolar waves with I = 2, m = 2 
and a stronger amplitude wave, with A = 0.001, cor- 
responding to a perturbation of 3% in the metric com- 
ponents. In this full 3D case, we do not use an octant 
of the spacetime, but rather set up a full grid with the 
origin in the center. Again, the grid has 120 3 points with 
Aa; = Ay — Az = 0.05. As we have a full grid, the outer 
boundary is now closer. In Fig. [l3| we see the evolution 
of the now stronger initial packet propagate outwards, 
as expected. In Fig. O, we again show the "waveform" 
measured directly by the function r(g xx — 1) at the outer 
boundary, which is allowing the wave to cleanly propa- 
gate off the grid. At late times boundary effects become 
visible. See Ref. for an excellent discussion of the 
problem of outgoing conditions in this scenario and pos- 
sible solutions with perturbative techniques. 



1.04 r 



1.02 



1.00 




0.98 L 



FIG. 13. We show the evolution of g xx for the I — 2, m = 2 
stronger amplitude quadrupolar wave packet along the 2-axis. 
The perturbation on the metric components is around 2.5% 
for this higher amplitude. Although this packet is fully 3D 
and can not be evolved using an octant of the spacetime, the 
metric component g xx is symmetric around the origin. 
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FIG. 14. We again show the time evolution of r(g xx — 1) at 
the outer boundary for the simulation shown in the previous 
figure. Again, although our simple copying boundary condi- 
tion, coupled with the MacCormack method, does a reason- 
ably good job of allowing the wave to propagate through the 
boundary, at late times boundary effects are evident. Note 
that the outer boundary, at z = ±3, is now closer to the 
origin. 

To better visualize the temporal evolution of this wave, 
in Fig. [l| we show the value of r(g xx — 1) along the 
z— axis evolving in time as a surface. We can see that 
the wave propagates cleanly away from the center and 
off the boundaries, as expected. 





FIG. 15. We show the time evolution of the "extraction" 
function r(g xx — 1) along the x line. The surface plot has time 
along the j/-axis. The r factor corrects for the 1/r fall-off, so 
we can see that the wave propagates from the center and off 
the boundaries. The previous figure corresponds to the y-axis 
(i.e., time) boundary of this plot. 
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The best way to visualize the full 3D nature of these 
waves and their propagation would be to show a movie, 
which obviously we can not do in printed form. In Fig. [l6| 
we show four snapshots of such a movie, showing two iso- 
surface values of the metric component g rr , constructed 
from the cartesian metric. 
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FIG. 17. We show the time evolution of the global conver- 
gence rate (computed in the £2 norm) of the hamiltonian con- 
straint, g xx and the lapse function. All quantities converge at 
second order. In particular, the hamiltonian constraint con- 
verges also at second order, although it does not converge to 
zero, since it is only satisfied to linear order. 



FIG. 16. We show four time snapshots of the evolution 
of the packet presented in the last figures. Two isosurface 
values of the spherical metric component g rr , reconstructed 
from the evolved cartesian components, are shown at times 
t = (a), t = 1 (b), t = 2 (c) and t = 3 (d). The dark 
and light colored isosurfaces correspond to the values 0.9997 
and 1.0003 respectively. They oscillate around the center and 
propagate outwards. 

All the wave tests presented in this section converge 
as expected. In Fig. 17] we show the time evolution of 
the convergence rate a obtained using the £2 (i.e., RMS) 
norm over the entire grid. We measure convergence for 
the hamiltonian constraint, the metric component g XXl 
and the lapse, and note that all converge at or above 
second order, as expected. We do not need to resort 
to the linear solution to measure convergence. In this 
special case, we do not measure the convergence of the 
constraints to zero, as the initial data is linear and only 
satisfies the constraints to linear order. Thus, we measure 
the convergence of the hamiltonian as we would any other 
quantity, using three different resolutions to create the 
convergence exponent. 



VI. BLACK HOLE TESTS 

Black hole spacetimes are currently one of the major 
motivations for developing 3D numerical relativity. The 
waveforms emitted by inspiraling colliding black holes 
are expected to be one of the most likely candidates for 
early detection by laser interferometers ||[|, and hence 
are urgently in need of general 3D simulations. Thus, 
black hole spacetimes are important tests of our code, 
and we will follow the work of Ref. |l9| in these code tests 
of Schwarzschild black holes. More dynamic black hole 
studies, including the simulations of 3D excitation and 
ringdown of the quasinormal modes of distorted black 
holes f^ , |7^J^ , |96f , and of black hole collisions jig] , are 
in progress and will reported and compared against pub- 
lished results in a future paper in this series. 

Black hole spacetimes are in many ways similar to 
other spacetimes. An initial metric evolves with some 
slicing conditions, and the constraints should converge 
as in any spacetime. However, special difficulties are en- 
countered due to the presence of singularities. Thus, as 
discussed in the introduction, present Cauchy evolutions 
of general 3D black hole spacetimes do not allow a 3D 
code to run forever, as they can when propagating dis- 
turbances in flat space or low amplitude waves. At some 
point, a time slice may hit a singularity and crash, or 
stretch the grid so much that the simulation will no longer 
be able to continue. At this point, we will see "blow ups" 
on our grid, convergence will fail (starting, usually, at the 
lower resolution grids), and we will have to stop our code. 
Thus evolving black holes for many tens of M, where M 
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is the ADM mass, with a demonstration of convergence 
is still "state of the art" in numerical relativity. 

In this section, we test Cactus using a single black 
hole with the Einstein-Rosen bridge topology with an 
isotropic radial coordinate r. That is, the spatial line 
element takes the form 



ds 2 = V 4 (dr 2 +r 2 dn 2 ) 



with 



# = 1 



M 
2r' 



(59) 



(60) 



We satisfy the constraint equations with this metric and 



initial Kij = 0. For more detail, see Ref. 19 1. For all the 



work which follows, we choose M = 1. 

This data is isometric in inversion through the sphere, 
or throat, located at r = 4p The singularity at r = is 
also related to the remapping of a second universe on the 
other side of the bridge to the origin in our flat space. 
However, rather than evolve the Einstein-Rosen bridge 
black hole spacetime with the natural S 2 x R topology 
(as used in axisymmetric simulations such as |Mp8fl ), 
we evolve it on an i? 3 manifold which contains a point 
where the conformal factor is infinite. This was one of 
the techniques used in Ref. Jn|, and has recently been 
generalized to generate full 3D, binary black hole data 
with spin and momenta (9^]. 

As in Ref. (T^j2^], we handle the infinity in the con- 
formal factor numerically using two tricks. First, we do 
not place a grid point at r — 0, but rather we stagger 
the origin, with grid points at Ax/2 and —Ax/2. Sec- 
ondly, we exploit knowledge of the conformal factor and 
its derivatives in our finite differencing. This allows us 
to factor out the infinity from the evolved quantities as 
known derivatives in the source terms, and evolve fields 
which are unity everywhere. This approach to comput- 
ing "conformal derivatives" is quite general, and can be 
used with a numerically generated initial data set as well. 
Note that this c onfor mal rescaling of the equations, as 
discussed in Sec. II D and Appendix A. is different from 
the conformal rescaling done in typical ADM codes, in- 
cluding the Cactus ADM thorn, where the Ricci tensor 
is formed directly with conformal derivatives of the sys- 
tem. For our first order system, we do not form the Ricci 
tensor, and therefore we must treat the conformal rescal- 
ing differently in order to preserve a first order system, 
and still allow only conformal variables to appear in the 
fluxes. 

Here we consider various slicings of a single black hole 
spacetime. We do not discuss or demonstrate multiple 
black hole or distorted black hole spacetimes here, since 
we wish only to show code tests at this time. However, 
preliminary tests show that the results presented here 
carry over into more dynamical black hole spacetimes. 
This is a major and active research area in 3D numeri- 
cal relativity in which we arc presently engaged. In the 
final part of this section, we also perform tests of the 



Schwarzschild black hole system with the ADM equa- 
tions in Cactus, and compare with the results from the 
BM formulation. All simulations in this section are done 
with Ax — Ay — Az, nx = ny = nz, and with the 
conformal rescaling of the BM system, or conformal dif- 
ferencing in the ADM system. For the BM system, all 
simulations were performed with the Einstein system. 



A. Geodesic Slicing 

A black hole spacetime evolved with geodesic slicing 
(a = 1, (3 l = 0) can only be evolved until points initially 
on the throat hit the singularity unless points are excised 
from the grid, as shown in Refs. |l^j2^]. At that point 
any code evolving this system will crash. We know that 
observers initially at rest in the Schwarzschild spacetime 
that this crash must come at t = irM. The crash will 
appear as an infinity or undefined value at a point on the 
numerical grid. 

Despite this critical failing, the geodesically sliced 
Schwarzschild spacetime is useful as an analytic solution 



for the three-metric exists, the Novikov solution [100 



This solution expresses the metric in terms of cyclic in- 
fall times for initially non-moving observers. Expressions 
for these solutions in isotropic radius are given in [ p0[ , 
although the final term is missing a squ are r oot, thus for 
completeness we give expressions here 
slightly different notation than Ref. 
radius, r„ is the areal radius, r„ 



101]. We use a 



r is our isotropic 
is the maximum (areal) 
radius for an observer during the cyclic infall (and is 
therefore the initial areal radius, so r a = r max at r = 0), 
r is the proper time of an observer (and therefore the 
coordinate grid time, as a = 1), and g rr and ggg are the 
conformal isotropic metric components. The relevant ex- 
pressions for an M — 1 black hole are 



<(r) = 



(1 + 2rf 



4r 



' a > ' max } 



L± 1 



dr n 



«(¥)' 

3 _ r a 

3 / ^max 




(61a) 



(61b) 



9rr 

and 
gee - 



2 V r a 

dr a 
dr m£ 



1 arccos 




(61e) 



To const ruct the metric, we must numerically invert rela- 
tion Eq. (|61bD to find r(r, r max ). Simple bisection solves 
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this problem. Aside from this minor complication, con- 
structing the solution is straightforward. 

We present here two demonstrations that our code is 
in fact creating the correct solution for the geodesically 
sliced black hole spacetime. In Fig. [if] we show the dif- 
ference between the g rr produced by the code (which is 
constructed from the full evolved cartesian three metric) 
and the analytic expression in Eq. (^3l]). We extract the 
data along a diagonal line. We show the difference at 
three different resolutions, adjusting the lower resolution 
differences by factors of 1/4 and 1/16, respectively. We 
note that the points (shown as crosses, diamonds, and 
triangles) are, for all practical purposes, identical in this 
figure, strongly indicating second order convergence at 
every point on the grid. 
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FIG. 18. We show the difference of the radial met- 
ric between the analytic Novikov solution and the full 
three-dimensional numerical evolution. Data is extracted 
along a diagonal line. We define E as the difference between 
the analytic solution and the numerical solution. We show 
E/16 for Aa; = 0.2, E/4 for Ax = 0.1 and E for Ax = 0.05. 
We note that the data points are practically identical, showing 
second order convergence. 

In Fig. [l9|, we show similar plots for the hamiltonian 
constraint. We show the lower resolution constraints di- 
vided by 4 and 16 respectively. Once again, the fact that 
these lines are visually coincident (although not com- 
pletely identical) strongly demonstrates that our code is 
converging at second order. Note that the error is largest 
near the throat, located at r = 0.5M, which is well in- 
side the horizon at late times as it rushes towards the 
singularity. 
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FIG. 19. We show the hamiltonian constraint, H, for the 
geodesically sliced black hole at three different resolutions. 
We show H/16 for Aa; = 0.2, H/4 for Ax = 0.1 and H for 
Ax — 0.05. We note that the lines are identical, indicating 
second order convergence. 

From these two diagrams we can calculate four val- 
ues of the convergence exponent everywhere, since we 
have two quantities to measure against exact solutions, 
at three resolutions. Doing this analysis gives a conver- 
gence exponent between about 1.8 and 2.1 (oscillatory 
in time), once again demonstrating second order conver- 
gence to an exact solution. 

The black hole spacetime also provides a strong test 
of a code's ability to preserve the appropriate spherical 
and rotational symmetries inherent in the initial data set. 
Especially near the singularity, there is a rapid growth of 
strong gradients surrounding a black hole, which must 
be computed in the separate cartesian metric functions 
on a cartesian grid. These individual functions do not 
exhibit the underlying symmetries of the black hole, so it 
can be difficult to model spherical or axisymmetric phe- 
nomena without introducing spurious effects due to reso- 
lution and coordinate geometry. As noted, the MacCor- 
mack method exactly obeys rotational symmetries, with 
g xx along an x-line through the origin and g yy along a 
y-line through the origin being the same to machine pre- 
cision, but has no such property for spherical symmetry. 
Therefore we must test whether spherical symmetry is 
preserved. We can make this test visually, by display- 
ing all the points on our grid versus their radial position 
for some spherically symmetric quantity, such as g rr in 
a "scatter" plot. In Fig. ^o] we do exactly this for the 
high resolution geodesically sliced black hole above. The 
"width" of the scatter plot at late times indicates that 
deviations from sphericity are becoming larger as the so- 
lution evolves towards the singularity. 
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FIG. 20. We use a scatter plot to show g rr vs. r for 
all points in a region of the cartesian grid for a geodesically 
sliced black hole. This scatter plot allows one to see how well 
spherical symmetry is maintained by eye. From this plot, it 
is clear that deviations from sphericity occur near the peak 
at late times, and are fairly small. This figure was generated 
with Ax = 0.05M, and the slices are shown at t — 0.9M, 
1.8M, and 2.7M. 

We repeat these tests using our Lax-Wendroff direc- 
tionally split update method. In Fig. |2l] we show a scat- 
ter plot of the conformal radial metric function in the 
neighborhood of its peak at r = M/2. We notice that 
spherical symmetry is obeyed very well despite the fact 
that this method is a manifestly cartesian method. 



FIG. 21. We use a scatter plot to show g rr vs r for all 
points in a region of our cartesian grid for a geodesically 
sliced black hole. This simulation uses the directionally 
split Lax-Wendroff solver. This figure was generated with 
Ax = 0.05A4" and slices are shown at t — 0.6M, 1.5M, and 
2AM. We note that even though the method is explicitly 
split in cartesian directions it maintains excellent spherical 
symmetry. 

In Fig. [2^ we show the solution is indeed converging 
at second order. We measure the convergence order of 
the hamiltonian constraint and find it is converging at or 
above a — 2 during the entire evolution. 



io"; 
io" : 



10" 
io" ; 



10"' 





A H (Ax 0.1) 




H/4 (Ax 0.2) | 






! t=0.8M W 




















A 




! t=2.4M 







10" 

10- 4 

10" 5 
10" 6 

10- 7 

01234567 
x/M 

FIG. 22. We show the convergence of the hamiltonian con- 
straint for a geodesically sliced black hole with the flux evolu- 
tion solved with the directionally split Lax-Wendroff method. 
We note that convergence is excellent, and points away from 
the boundary are visually coincident, demonstrating second 
order convergence of the constraints towards zero away from 
the boundaries. 



B. Algebraic Slicings 
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Algebraic slicing conditions have been used for three 
dimensional black hole evolutions in the past with a rela- 
tively high degree of success, as shown in Refs. |l9| , |2"3| , p4j . 
Such slicings typically use Eq. (|6|) to provide a condition 
on the lapse. Here we examine the use of such slicings 
in 3D black hole spacetimes in the BM formulation. We 
note that such slicings also have been shown under cer- 
tain conditions to develop coordinate pathologies p5[ , 
but we will not investigate those issues here. The main 
purpose of this section is to compare results of Cactus 
with previously published results on Schwarzschild black 
hole evolutions. 

The simulations in Refs. |l^,|23|,^4j used both a dif- 
fusion term added to the lapse evolution equation to 
achieve stability, and an enforced isometry condition, 
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mapping the highly resolved region exterior to the throat 
into the poorly resolved region interior to the throat. As 
detailed in Ref. jL9| , explicit enforcement of this isometry 
was very important in obtaining accurate long time evo- 
lutions of the system, as it allows one to avoid numerical 
evolution in the coarsely resolved region near the singu- 
larity: one simply maps the accurately computed exterior 
into this region before proceeding to the next time step. 
Although the algebraic slicing conditions studied actu- 
ally do obey the isometry operation, and will attempt to 
preserve it numerically during an evolution, without an 
explicit isometry operator in the code, large errors will 
develop inside the throat, causing a code to crash. 

An isometry condition could be applied within Cactus, 
but with the BM system this leads to an additional dif- 
ficulty in that the isometry conditions on the Dijk, A4 
and Vi variables is non-trivial, since these are not ten- 
sor quantities. Due to this complication, and due to the 
promise of alternative techniques such as apparent hori- 
zon boundary conditions which do not require isometry 
conditions, we have currently chosen not to implement 
an isometry condition in Cactus. Under these conditions 
it is difficult to achieve the same accuracy and long run- 
time that were available to an isometry based code, when 
algebraic slicings are used. We stress that this not a lim- 
itation of the code or the formulation of the equations, 
but merely a sensitivity of such slicings in black hole sim- 
ulations without an explicit isometry operator. Similar 
results are obtained, for exa mple, w ith the "G" code used 
to generate results in Refs. [ ^9|j2^j2^] . Furthermore, we 
will see in the next section that maximal slicing, which as 
shown in Ref. |fl9|| does not require the isometry operator, 
works very well in Cactus. 

With those remarks in mind, in Fig. ^3] we show the 
lapse profile in a scatter plot at t = 2.1M and 3.5M 
for a "1 + log" (/ = 1/a in Eq. @) sliced black hole. 
Clearly the failure of the lapse to collapse in the center, 
combined with poor resolution of the consequent gradi- 
ent in the lapse, will not allow an accurate evolution in 
this region. Convergence testing this solution at three 
resolutions, Ax = QAM, Ax = Q.2M, and Ax = QAM, 
in this case on a full grid rather than an octant grid, we 
see that the simulation on the medium resolution grid 
crashes first, at t = 4.2M. We note that the lapse is col- 
lapsing most quickly at r — M/2, indicating that our sys- 
tem is trying to preserve the underlying isometry present 
in the initial data, as it should. 
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FIG. 23. We use a scatter plot to show the collapse of the 
lapse at t = 2.1M and t = 3.5M for a "1+log" sliced black 
hole. We note that with this dynamical local slicing, spher- 
ical symmetry is not preserved to such a high degree as in 
the geodesically sliced case, especially near the point of large 
gradient in the lapse. This inaccuracy in the dipping of the 
lapse will cause the code to crash shortly after this plot. This 
data was produced with Aa; = 0.1M. 

We can understand the nature of the algebraically 
sliced spacetime without isometry or diffusion by study- 
ing its convergence properties. In Fig. |24| we show the 
convergence of the hamiltonian constraint towards zero 
at three different resolutions. Several important points 
in this figure should be noted. Firstly, the very low 
resolution simulation (Ax = QAM, almost the radius 
of the throat) does not converge from the second slice, 
t = 2.5M. As we shall see below, this is because this very 
low resolution simulation "misses" the isometry, and the 
lapse collapses, leading to the longest time evolution of 
the three simulations. We note, however, that at the first 
two displayed times, the medium and high resolution sim- 
ulations are converging appropriately, as indicated by the 
almost coincident lines in Fig. ^J. As the medium res- 
olution grid nears its crash time at t = 4.2M, however, 
there is no strong evidence of second order convergence 
in the system near the hole. 
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FIG. 24. We show the convergence of the hamiltonian con- 
straint for the "1 + log" sliced spacetime. We show 1, 4, and 
16 times the hamiltonian constraint at Aa; = 0.4M, 0.2M, 
and 0.1M respectively. We note that the low resolution grid 
does not converge, and the entire system fails to converge at 
late times. The medium resolution grid simulation will crash 
at t = 4.2M. 

Returning to the mystery of the lowest resolution grid, 
we show in Fig. ^5] the evolution of the lapse on the low 
(Ax — 0.4A/) and high (Aa; = 0.1M) resolution grids. 
On the lowest resolution grid the system simply has too 
few points to obey the isometry, and the lapse collapses 
uniformly around r = 0, allowing a long time evolution. 
The highest resolution grid clearly attempts to obey the 
isometry, but is destined to fail, due to the small number 
of points covering the region in r < M/2. Thus the 
two evolutions do not converge towards the same slicing 
of the spacetime. This is clearly a dangerous feature 
running a simulation with very poor resolution: it can 
produce a solution which misses features, but still creates 
a reasonable looking (and, in this case, longer running) 
solution than a higher resolution run. 

To summarize this subsection, algebraic slicings are 
convenient and inexpensive singularity avoiding slicing 
conditions. We have shown that in Cactus they behave 
as expected, and converge at second order as they should. 
But they must be used with care. As already shown in 
Rcf. p9[ , in a spacetime containing a singularity, they 
can still be very useful if an isometry operator is used to 
avoid evolving the region near the singularity. Without 
it, evolution in this region inside the black hole throat is 
almost impossible in 3D. These slicings should still find 
use in a number of other circumstances, including use on 
black holes if precautions are taken near a singularity, if 
it exists on the grid. 
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FIG. 25. We note that on a low enough resolution grid, the 
isometry condition inherent in the Schwarzschild spacetime, 
which maintains a = 1 at the center of the grid, is "missed" 
so the lapse will collapse everywhere in the center of the grid. 
Ironically, this under-resolution will allow the low resolution 
simulation to run considerably longer than the high resolution 
simulation, but the solution will not converge, of course. Lines 
are shown with Aa; = 0.1M (solid line) and QAM (dashed 
line) at t = 0, 1, 2, 3, and AM. 



C. Maximal Slicing 

Maximal slicing has long been a favorite slicing condi- 
tion for numerical relativity. Alas, the maximal slicing 
condition, Eq.(p8|), is an elliptic condition for the lapse, 
which brings with it both a breaking of the hyperbolic 
system for the lapse and its derivatives and a very large 
degree of computational complexity. 

Solving three-dimensional elliptic equations is far more 
difficult than solving three-dimensional hyperbolic ones, 
using much more memory, and taking much more time. 
For the work in this paper, we use an elliptic solver based 
on the freely available PETSc software [ 102| , which uses 
Krylov subspace based matrix methods, such as conju- 
gate gradient, to solve the elliptic conditions which are 
rewritten as a matrix equation after being cast in finite 
difference form. The Cactus code has several additional 
elliptic solvers with various degrees of efficiency and func- 
tionality, including several relaxation based solvers, and 
a parallel multigrid solver developed by B. Briigmann, 
based on the solver used in Refs. fl9^j46| . 

There are various boundary conditions we can apply 
to the lapse at the outer boundary when using maximal 
slicing. For example, we can allow the boundary value of 
our lapse to change in time, applying the same boundary 
condition to both the lapse and its derivatives that we 
apply to all other fields. The "flat" boundary condition 
used here has the effect of copying the lapse from one 
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point in the interior to the exterior after the maximal 
equation solve, which causes the lapse to collapse slowly 
at the outer boundary. We can also use the more tradi- 
tional approach, which keeps the boundary fixed at some 
initial value for the entire run. Experience has shown 
that the best approach in this case is to be initialize a to 
the (static) Schwarzschild value, and then call the max- 
imal solver, to create the initial lapse profile with the 
correct spherical outer boundary, as discussed in Jl9| . 
This has the added advantage of holding the lapse at the 
Schwarzschild value near the boundary, reducing evolu- 
tion of the metric there for some time. 

As shown in Ref. unlike with algebraic slicings, 
one can handle the region inside the throat of the black 
hole simply by ignoring it. The elliptic maximal lapse was 
found to collapse rapidly throughout this troublesome re- 
gion, quickly halting the evolution there. Hence no spe- 
cial precautions, and no isometry operator, are needed to 
handle this region. We will see the same behavior in Cac- 
tus below. Although maximal slicing could be enforced 
with an isometry condition, as in axisymmetric simula- 
tions |^3,103,37 104 1, it is not necessary to do so, and we 
shall not do so here. 

We show here that Cactus runs and converges using 
maximal slicing by evolving a single black hole on a 100 3 
and a 50 3 sized computational grid to 15M. We look for 
convergence in the constraints, which should converge to 
zero, and also in trK, which the maximal slicing condition 
should force to zero. 

First we show the behavior of the solution. In Figs. |26| 
and [27] we show the "collapse of the lapse" and "peak 
in g rr " which are familiar from maximally sliced black 
hole evolutions in numerical relativity. Due to the sin- 
gularity avoiding coordinate slicing, we see the lapse col- 
lapses towards zero at the center of our grid, which halts 
evolution, and leads to large proper distances between 
coordinate points as the exterior evolves, creating large 
gradients in the metric. We also note that the lapse does 
collapse at the outer boundary in this simulation, and 
also that at late times, outer boundary effects are notice- 
able in g rr /ip A . 




FIG. 26. We show the collapse of the lapse along the :r-line 
for a maximally sliced black hole. We note the traditional col- 
lapse in the center. We also note that our outer boundary is 
not held static in this case, and thus the lapse collapses there. 
This collapse allows evolution with the outer boundary placed 
nearer the hole than in the static boundary case. This simu- 
lation has Ax = 0.1M and the lapse is shown every t — Q.8M 
from t = 0M to t = 14.4M. 
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FIG. 27. We show the growth of the conformal 3-metric 
g rr /tjj 4 along the a;-axis in the maximal slicing case. This fig- 
ure is the now infamous "grid stretching" figure, and demon- 
strates the problem which plagues all black hole simulations 
with singularity avoiding slicing without apparent horizon 
boundary conditions, namely the explosive growth of the ra- 
dial metric function. Late time outer boundary problems are 
also evident in this plot. This simulation has Ax = 0.1M 
and the metric is shown every t = 0.8M from t = 0M to 
t = 14.4M. 
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In Fig. g8| we demonstrate that the system maintains 
spherical symmetry using the initial lapse of one and al- 
lowing the lapse to change at the outer boundary. We 
note that even on this log scale, and at very small values 
of the lapse (a — > 10 -4 ) the system maintains excellent 
spherical symmetry 
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FIG. 28. We use a scatter plot to show the maintenance of 
spherical symmetry in the lapse in a maximally sliced black 
hole. We show the lapse on a log plot, and note that the 
collapse to very small lapse maintains spherical symmetry to 
a very high degree. Slices are shown at t — 3.6M, 7.2M, 
10.8M, and 14AM The resolution used is Ax = 0.1M. 

We emphasize that the growth in g rr is not something 
we simulate directly. We do not evolve g rr , but rather, 
we evolve cartesian metric functions. These functions 
must not only display the growth in the radial metric 
function, but must also contain the decreasing behavior 
of the angular metric functions. That is g xx must behave 
like g rr along the x— line, but also like gee along the y— 
and z— lines. This leads to an even larger dynamic range 
in our cartesian metric functions than in the radial or an- 
gular metric functions alone. In Fig. we demonstrate 
this by showing a slice in the x — y plane of g xx ji\)^ at 
t = XAAM for the high resolution (Ax = 0.1M) simu- 
lation considered above. It is clear from the figure that 
the function is growing along the x— axis and dropping 
along the y— axis, as expected. 




FIG. 29. We show the behavior of the cartesian metric 
functions by showing a slice in the x — y plane of g xx /ip 4 
at t = 14 AM for the high resolution (Ax = O.lJVf) simulation 
considered above. Note that this function behaves like g rr 
along the x— axis, and ggg along the y— axis. 

One obvious question to ask of our simulation is 
whether or not our maximal slices are actually maximal, 
in that they keep the trK zero. This condition will be 
violated by our numerical simulation, but we can check 
whether the trK converges towards zero. In Fig. ^ we 
do exactly this. We show the tiK/4 on the 50 3 Aa; = 0.2 
grid and trK on the 100 3 Aa; = 0.1 grid. If the trK 
converges to zero at second order we expect these lines 
to be identical. From Fig. [?0] we can clearly see that our 
evolution converges to a maximal slice at second order. 
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FIG. 30. We show the trK at two different resolutions for 
the maximally sliced case. The dotted line is trK/ 4 calcu- 
lated with Ax — 0.2A/. The solid line is the trK calcu- 
lated with Arr = 0.1M. Since the trK = should remain 
constant in maximal slicing, plotting these two quantities to- 
gether demonstrates the second order convergence of our so- 
lution to the maximally sliced spacetime. 

Similarly we can confirm that we are creating a solu- 
tion to the Einstein equations in our maximally sliced 
spacetime by testing if the Hamiltonian constraint con- 
verges to zero. In Fig. |3l] we show one-quarter of the 
constraint on the 50 3 Ax = 0.2 grid and the constraint 
on the 100 3 Aa; = 0.1 grid. Once again, we see the lines 
are close to identical visually, strongly indicating second 
order convergence. 
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FIG. 31. We show second order convergence of the hamil- 
tonian constraint in the maximally sliced black hole. We show 
H at Ax = 0.2M and H/4 at Ax = 0AM. The fact that the 
points are visually coincident strongly indicates the system is 
converging towards a solution of the Einstein equations. 

Figs. ||c] and [^demonstrate fairly conclusively that in 
the regions where our error is larger, our code is converg- 
ing at second order. However, we note that our simple 
boundary conditions do lead to small, but non second or- 
der con verg ent, errors at low levels which are not visible 
in Figs. 30 and |3l]. We can see these by plotting the trK 
and hamiltonian at late time (here t — 14AM, the final 
time in Figs. p0| and pl|) u sing a logarithmic y axis, which 
we do in Fig. 133. Fig. |32| artificially inflates the non sec- 
ond order convergent features due to the boundary, but 
it is instructive nonetheless. Since these features are on a 
very low level (several orders of magnitude smaller than 
the dominant error) they have no real adverse affect on 
our solution at this time. 




lO" 1 



,-2 



10 
1(T 



10" 
10 



10 



,-5 



r %. 










r trK (Ax 0.1y*******> 


^^^^^ ! 






I A A trK/4 (Ax 0.2) 


* •> A i 

A 1 



10 



x/M 



FIG. 32. We show the convergence of trK and the hamil- 
tonian constraint at t = 14AM for the single black hole case 
considered in Fig. |3o| and |3l| Here we use a logarithmic y 
axis, which emphasizes that, at a very low level, the bound- 
ary introduces non second order convergent features to the 
system. Since these effects are several orders of magnitude 
below the dominant errors, they do not have an adverse effect 
on our solution (the area with large error is what crashes our 
code). However, it is clear that our boundaries lead to small, 
but non second order convergent, effects entering our grid in 
the less dynamic (spatial) regions. 

For the sake of completeness, we also show the fail- 
ing of convergence of these quantities when our grid is 
too poorly resolved. In Fig. |33j we show the constraints 
for a 25 , Ax = 0.4 run. Since only two points cover the 
entire initial horizon at this resolution, we cannot reason- 
ably expect a converged answer, and we see that, even 
though the higher resolution simulations are converging 
at second order, the low resolution simulation has a worse 
convergence property. At late times, this effect is mostly 
due to the lower resolution leading to an earlier crash 
time on the lowest resolution grid. 
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FIG. 33. We repeat the display in Fig. |3lj, adding a third 
very low resolution simulation with Aa; = AM. We see that 
the low resolution solution is not converging for t > 12M, and 
it will crash soon after the final slice shown here. 

We repeat these tests with the outer boundary on the 
lapse held fixed, which is the condition used previously 
in Ref. |pj| . For this case, we must set the boundary 
farther away than in the non-fixed case, setting it here 
initially at t — 15M, rather than t = 10M. Additionally, 
we must use the maximal slicing solver to generate an 
initial lapse which has the isotropic Schwarzschild form 



2r — M 
2r + M 



(62) 



at the outer boundary. This leads to an initial lapse 
other than one everywhere, as discussed in Ref. p9| . By 
holding the lapse fixed, we avoid the collapse of the lapse 
near the boundaries, and therefore evolve for a somewhat 
longer proper time in the outer region. (We note that 
with a boundary very far away, as could be provided by 
some form of adapted mesh structure, the two conditions 
would be equivalent). Despite this difference, when we 
evolve the maximally sliced black hole system with the 
two boundary conditions, we see quantitatively the same 
behavior in the metric functions. 

In Fig. [m| we show the lapse along the x line up to 
t = 14AM, and note it remains fixed at its initial outer 
boundary value, as expected. Comparing with Fig. |2^, 
we can clearly see that this stops the lapse from collapsing 
over such a wide portion of the grid, with the lapse at 
x = 10M being around 0.9 in Fig. |34|, and closer to 0.7 
in Fig. p6|. 




FIG. 34. We show the collapse of the lapse in the maxi- 
mally sliced spacetime with the outer boundary held fixed. 
Note the initial lapse is not one, as it was in Fig. but 
rather was the result of applying our maximal slicing solver 
to a lapse which obeys the Schwarzschild lapse outer bound- 
ary conditions. The resolution used here was Ax — 0.15M. 
Slices are shown every 1.2M between t — OM and t — 14AM. 

In Fig. |3| we show the convergence of the tvK and 
hamiltonian constraint to zero at t = 12M, and note we 
still achieve second order convergence visually. 
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FIG. 35. We show the convergence of the trK and the 
hamiltonian constraint towards zero at t — 12M for the max- 
imally sliced spacetime with a fixed outer boundary. We 
note that, due to the lower resolution used in this simulation 
(Aa; = 0.15 rather than 0.1) convergence of the constraints 
starts to fail near the peak earlier. 
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D. Comparison with ADM Code 

In the above sections, we have demonstrated that the 
BM formalism can generate convergent black hole space- 
times. In this section we confirm that the Cactus ADM 
integrator is also second order convergent on black hole 
spacetimes by repeating several of the above tests, and 
compare ADM results with results from the BM formu- 
lation. 

In general, we find that the BM and ADM systems 
generate comparable results, although, as shown below, 
the ADM system we have implemented will generally run 
some time longer than the BM system in maximally sliced 
black hole spacetimes, with large errors appearing first 
in the BM system. The grid stretching problems ulti- 
mately ruins both calculations. We emphasize that this 
is in not a shortcoming of the BM system; treating the 
system with advanced methods as described in, for in- 
stance, Ref. will allow the BM system (in one di- 
mension) to evolve for significantly longer times than the 
ADM system, showing a real advantage in using the first 
order system when combined with advanced numerical 
techniques. Rather, this demonstrates that when evolv- 
ing mathematically equivalent systems of equations, on 
problems such as these black holes that have large gradi- 
ents, without using numerical methods designed to han- 
dle such features, both will fail when gradients become 
too steep to resolve. The details of how the calculation 
fails can depend on many factors. Thus for black holes, 
the present numerical methods applied to to BM system, 
which only moderately exploit the first order nature (in 
this case, "flat" boundaries, the Strang split, and the 
true MacCormack method) are not guaranteed to gener- 
ate markedly better or longer numerical evolutions than, 
say, the ADM system with leapfrog. 

We return first to the geodesically sliced black hole. 
In Fig. ||(| we repeat the test from Fig. [ll| by compar- 
ing g rr hv with the analytic Novikov solution, Eq.(|6T]). 
We show the difference at a high resolution (Ax = 0.1) 
and 1/4 the difference at a lower resolution, (Aa; = 0.2). 
The points are visually coincident indicating second or- 
der convergence on our grid, which we see in general. 
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FIG. 36. We show that the ADM integrator in the Cactus 
code converges to second order against the analytic Novikov 
solution, repeating the test presented with the BM integrator 
in Fig. hq. We show the difference between the analytic solu- 
tion and the computed solution at Aa; = 0.1 and one quarter 
the error at Ax — 0.2. The fact that the points are visually 
coincident demonstrates second order convergence, which we 
see on our entire grid. 

The most interesting comparison is the maximally 
sliced black hole. Studies of the three-dimensional maxi- 
mally sliced black hole with the ADM system have been 
undertaken in great detail in Ref. p9|, so we only treat 
them briefly here, using the cactus ADM integrator. 

In Fig. [37| we show that the BM and ADM system 
give qualitatively the same behavior at a fixed resolution 
(since both systems converge at second order, in the limit 
of infinite resolution they would give identical results). 
We show the metric function g rr /ip 4 every 3M between 
0M and 15M with a resolution Ax = 0.15 and a 100 3 
grid. We see clearly that the behavior is the same in 
both cases, although at late times, the two solutions are 
noticeably different near the peak in g rr . 
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FIG. 37. We compare g rr /ip 4 for a maximally sliced single 
black hole spacetime evolved with Aa; = 0.15 on a 100 3 grid 
with the ADM and BM integrators. In the ADM system 
all values are held fixed at the outer boundary, while in the 
BM system, only the lapse and its derivatives are held static, 
corresponding to the run in Fig. We show data every 3M 
between and 15M along the a;— line. We note that both 
systems exhibit qualitatively the same behavior. 

Even though the two solutions in Fig. [37] are different, 
both solutions are converging to second order, as shown 
in Fi gs. B8| and 39| . In these figures, we repeat the tests of 
Figs. |30| and pi by convergence testing both the tiK and 
the hamiltonian constraint against zero. We see converge 
close to second order visually in the figures, and every- 
where on the grid when measured globally. 
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FIG. 38. We repeat the convergence test in Fig.|3l] with the 
ADM integrator. We use parameters Aa; = 0.15 on a 100 3 
grid and measure H at Aa; = 0.15 and H/4 at Aa; = 0.3. We 
see almost second order convergence visually, and measure a 
convergence exponent around two over our entire grid. We 
note the the convergence order drops away from two as we 
approach the end of our simulation. 
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FIG. 39. We repeat the convergence test in Fig.^ with the 
ADM integrator. We use parameters Aa; = 0.15 on a 100 3 grid 
and measure trK at Aa; = 0.15 and trA'/4 at Ax = 0.3. We 
see second order convergence visually, and measure a conver- 
gence exponent around two over our entire grid. 

We finally directly compare the BM and ADM evolu- 
tions of the maximally sliced black hole spacetime with 
parameter Aa; = 0.15 on a 100 3 grid, with all fields held 
fixed at the outer boundary in the ADM system, and 
the lapse held fixed with other fields having the "copy" 
boundary conditions in the BM system. We calculate 
the hamiltonian constraint using Eq. ( p2| ) in both the 
ADM and BM system, constructing the BM -Dyfc and 
Vi variables from the ADM system with centered finite 
differences. 

In Fig. [lO] we can see that for a large part of the run 
time, the hamiltonian constraint, although different, is of 
the same order, around 0.1. However as t — > 15M, the 
hamiltonian constraint for the BM system around the 
peak in g rr drops to a larger (absolute) value than the 
ADM system. This dropping continues, causing the BM 
system to crash about 4 — 5M before the ADM system 
with the parameters chosen here. 
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FIG. 40. We compare the hamiltonian constraint in the 
BM and ADM systems for the Ax = 0.15 100 3 maximally 
sliced black hole simulation. The hamiltonian constraint is 
evaluated by Eqj22| in both cases, with the BM Dijk and Vi 
variables constructed from the ADM simulation at every time 
step. We note that the errors in the constraint are compara- 
ble, but at late times, the errors in the BM system are larger 
near the maximum of the grid stretching. In the simulation 
shown here, the ADM code will run around 4—5M longer than 
the BM simulation (with crash times around 16M and 20M 
at this resolution). We note that the constraints converge to 
zero in both cases. 



E. One-D AH Finder as a test of spherical 
spacetimes. 

Since the only black hole spacetimes we treat here are 
spherical, we can use spherical expressions for the loca- 
tion of the apparent horizon extracted along constant ra- 
dial lines of the spacetime. Here we choose diagonal lines. 
We assume the spherical metric has the line element 



dl 2 = ip A (g rr dr 2 + geedfl 2 ) , 
so the outgoing normal has the form 
1 



:(1,0,0). 



^ 2 y/9rr 

We can evaluate the expansion, 

D a s a + K ab s a s b 
2 
r 



(63) 



(64) 




trK = 
_ 2K eg 

^969 



(65) 



everywhere along this line. The point where the expan- 
sion crosses zero defines the apparent horizon. By mea- 
suring ip^r 2 ggg I 'AM 2 there, we can monitor the horizon 
area, which should be identically 1 using this normaliza- 
tion. 



In Fig. [y]we show the evolution of the apparent hori- 
zon area up to 15M in the maximally sliced cases dis- 
cussed above. We can make a crude estimate of how well 
our horizon is converging by measuring the convergence 
exponent for its radial location versus time. Although 
this measure is plagued by oscillations, we see that, on 
the whole, we have better than second order convergence. 
As well as a spherical AH finder, determined by finding 
the zero of Eq. (poj) , the Cactus code has a parallelized 
implementation of Gundlach's Pseudo-spectral apparent 
horizon finder [ |105| . Applications of this AH finder dur- 
ing dynamic evolutions will be discussed elsewhere. 
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FIG. 41. We show the area of the apparent horizon for 
the dx = 0.1M maximally sliced black hole. The apparent 
horizon is extracted along the diagonal line. 



VII. SUMMARY 

Hyperbolic formulations of Einstein's equations have 
been proposed by a number of groups as a promising tool 
for numerical relativity |^-]^, |c^ , [73| , |67| . These reformula- 
tions of Einstein's equations have shown great strength in 
ID tests @|. Early versions of the BM hyperbolic for- 
mulation pfwere developed into a full 3D code and tested 
on dynamically sliced flat space |5q] , leading further to 
the development of the "H" code which was applied to 
3D gravitational wave studies [|IJ. A 3D version of the 
Abrahams et al. hyperbolic formulation |64[] is also cur- 
rently under development |39| . But these 3D codes have 
seen only limited development and application. 

In this work we have performed the first systematic and 
detailed numerical exploration of a 3D hyperbolic formu- 
lation of Einstein's equations on a number of spacetimes 
of broad interest in physics and astronomy. We devel- 
oped and tested a full 3D numerical code, called Cactus, 
which implements the recent and more general BM hy- 
perbolic formulation of Einstein's equations pJ3] . With 
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this code, we showed on various dynamically sliced flat 
space, black hole, and gravitational wave spacetimes that 
this formulation allows for numerical treatment that is as 
stable and accurate as the traditional applications of the 
ADM formulation. 

The Cactus code has a modular structure allowing for 
different formulations of the Einstein equations, includ- 
ing the ADM system, different numerical methods, and 
many different initial data, gauge, and analysis routines. 
Cactus is developed on an advanced parallel computa- 
tional infrastructure, achieving over 66GFlops/sec on a 
512 node Cray T3E supercomputer ]70| , |7l|] . In this paper, 
within Cactus we compared Strang split MacCormack 
and Lax-Wendroff methods, applied to the Bona Masso 
system, against a leapfrog implementation of the ADM 
system, and also against previous results obtained from 
two completely independent 3D codes (the "G" code, 
based on the ADM formulation, and the "H" code, de- 
scribed above). The numerical methods used were de- 
scribed in detail. 

For the 3D black hole spacetimes, we studied (a) 
geodesically sliced black holes, and compared with the 
analytic solution of Novikov, (b) algebraic slicings, which 
have good singularity avoidance properties, (c) and max- 
imal slicing, which is has traditionally been a preferred 
choice for numerical black hole evolution. On all tests 
with both the ADM and BM formulation, the code per- 
formed well, reproducing previous published results on 
spherical black hole evolution in 3D . 

For 3D pure gravitational wave spacetimes, Cactus 
was tested on the evolution of linearized quadrupole and 
plane waves against previous results obtained with the 
"G" and "H" codes, again reproducing results of exten- 
sive studies published previously fl4l| . Cactus was also 
tested with dynamically sliced Minkowski spacetimes, 
where quantities such as Riemann invariants were shown 
to converge to zero. 

We also discussed the importance of convergence tests, 
and detailed a number of techniques we have developed 
to test convergence of the code. We showed that Cactus 
is rigorously second order convergent, and we emphasized 
that convergence tests are important techniques for di- 
agnosing code errors. 

We emphasize that although this paper shows many 
successful applications of a 3D hyperbolic formulation of 
Einstein's equations, we have focussed on applying stan- 
dard numerical methods for flux conservative systems, 
and on showing that they perform as well as standard 
methods applied to the ADM system. We have not yet 
exploited the kinds of advanced numerical methods that 
can be applied to the eigenfields of a hyperbolic system. 
Such numerical treatments are ultimately one of the ma- 
jor motivations for using hyperbolic systems in numerical 
relativity. The application of numerical methods specif- 
ically designed for hyperbolic systems (e.g. TVD meth- 
ods |5^^^]) has produced vast improvements in ID stud- 
ies of black holes, and their applications in 3D will are 
under development. Advanced inner (e.g. on a black hole 



horizon) and outer (e.g. at the edge of a numerical grid) 
boundary treatments may also be possible through the 
use of the eigenfields. The present Cactus code provides 
an advanced parallel tool for developing and testing such 
methods. 

This paper is the first in an anticipated long series 
with many collaborators. There are many directions in 
which research with this code is proceeding. We are cur- 
rently working on evolution of multiple black hole space- 
times, evolution of strong gravitational waves, 3D appar- 
ent horizon boundary conditions, self-gravitating scalar 
fields, advanced numerical treatments of the character- 
istic system, and full general relativistic hydrodynam- 
ics, among other projects. We expect that future papers 
will build on this one, continuing to show careful com- 
parisons with analytic solutions, demonstrating rigorous 
self-convergence, and discussing the effects of boundaries 
and numerical methods. 

We plan to make the code used for the all calculations 
in this paper publicly available at some point after the 
publication of this paper, together with the parameter 
files used to produce the figures and additional color fig- 
ures and movies that provide more details than it is pos- 
sible to show in printed form. All this information and 
instructions on uploading the code will be located at the 



web address http : / / cactus . aei-potsdam.mpg.de/ 
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APPENDIX A: THE CONFORMALLY 
RESCALED BM EQUATIONS 

Here we detail the equations modifications necessary to 
take into account a static conformal factor for the metric. 
We will evolve a "conformal" metric 9 ij related to the 
physical metric by the conformal factor ip: 



9ij = ip A 9ij ■ 



(66) 



We will keep the same formal definitions for the BM vari- 
ables: 



Vi = D ir r - D r 



(67) 
(68) 



where indices are raised with g^ and lowered with g^. 

With these definitions, the physical BM variables re- 
late to the conformal ones by 



Dkij = ip A (D kij + 2ipk9ij) 
V l = V l + 4^ . 



(69) 
(70) 



We also introduce the following notation for the deriva- 
tives of the conformal factor 



A = — 



(71) 
(72) 



The Christoffel symbols relate by 



t%=T\ j+ 26^ + 26^-2^ 9ij , (73) 



and the Ricci tensor by 

Rij = Rij — Yij — Y k g^ , 
where we define 



(74) 



(75) 



with covariant derivatives computed using r^ -. Thus, 

= 2g^ k ^ k + 2{i>ij - ^ r r r tJ - 3^j) . (76) 



We can then derive a modified set of fluxes (note that 
the flux for the extrinsic curvature does not change): 



Flgij = , 
F k a = , 

+1/2 6? (Aj+2Vj-Dj r r ) 



FtKij = -/3 k Kij + a [ D% - n/2 V k 9ij 



(77) 
(78) 
(79) 



+1/2 J? {Ai + 2Vi-D^)], 

F k D kt] = -(3 r {D rij + 2ip r9i j) + a (K^/^j 4 - Slj ) , (80) 

F k A k = -(3 r A r + aQ , (81) 

F k V t = -[3 k {V t + 4^0 + B k - B k . (82) 

The modified sources are: 

S -9ij =~2 a (K lj /ip i - Sij) 

+2p r {D ri j+2il }r g i j) , (83) 
S-a = -a 2 Q + af3 r A r , (84) 
S-Kij = 2{K ir B[ + K, jr B[ - KijB r r ) 

+a [ - (4) i?y + {-2Ki k K kj + trK K,j)/^ 4 
—T ri T r k j + 2D ik D r k + 2Dj k D ri + T kr T r ij 
-{2D k k -A r ){Dij + DjD (85) 
+Ai(Vj 1/2 Dj k ) + Aj(Vi 1/2 D k ) 
+Aj(Vi- 1/2 D ik k )-nV k D ki j] 

+n/4 agij [ -D k rs T k rs + D k r r D ks s - 2 V k A k 
+ (K rs K rs - (trK) 2 )/^ + 2a 2 G 00 ] 

-Y^ + 2A l ip J + 2Aji>i 

+ gi j [{n-l)Y k -2A k ^ k ] , 
S-D ki j = , (86) 
S-A k = , (87) 
S-V, = a 2 G° t + a/^[A r {K\ - tr K 51) 

+K r s (D ir s 2D r l) K\{D r s s - 2D J) 

-2^(3K r k -trK8l)} (88) 

+2{B{ - 81 trB) V r + 2{D r f - 5° D j jr )B r s 

+4Bi r ip r - AtrBipi . 

Finally, our algebraic slicing condition becomes 

Q = f(a)trK /^ . (89) 



[1] J. Masso and P. Walker, (1998), in preparation. 
[2] C. Bona and J. Masso, Phys. Rev. Lett. 68, 1097 (1992). 
[3] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 
Lett. 75, 600 (1995). 



33 



C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev 
D 56, 3405 (1997). 



Eanna E. Flanagan and S. A. Hughes, gr-qc/971012£ 
(1997). 



Eanna E. Flanagan and S. A. Hughes, gr-qc/970103£ 
(1997). 

A. M. Abrahams et ai, Physical Review Letters 80 



1812 (1998), kr-qc/970908l 

G. B. Cook et ai, (1997 ) |gr-qc/971107l- 
R. Gomez et ai, (1998), |g£qc/9801069|. 



G. Mathews and J. Wilson, Astrophysical Journal 482, 
929 (1997). 

M. Ruffert, M. Ram pp, and H.-T. Jan ka, Astron. As- 
trophys. , to appear, astro-ph/9611056 
M. Ruffert, H.-T. Janka, W. Keil, and G. Schafer, in 
Proc. of the 17th Texas Symposium on Relativistic As- 
trophysics (World Scientific, Singapore, 1995). 
K.-I. Oohara and T. Nakamura, in Relativistic Gravita- 
tion and Gravitational Radiation, edited by J.-P. Lasota 
and J. -A. Marck (Cambridge University Press, Cam- 
bridge, England, 1997). 

M. Choptuik, Phys. Rev. Lett. 70, 9 (1993). 

C. Gu ndlach, To app ear in Adv. Theor. Math. Phys 

(1998), lgr-qc/9712084 . 

A. Abrahams and C. Evans, Phys. Rev. Lett. 70, 2980 

(1993) . 

C. Gundlach, (1998), ^r-qc/9710066 . 

P. Anninos, J. Masso, E. Seidel, and W.-M. Suen, 

Physics World 9, 43 (1996). 

P. Anninos, K. Camarda, J. Masso, E. Seidel, W.-M. 
Suen, and J. Towns, Phys. Rev. D 52, 2059 (1995). 

B. Briigmann, Phys. Rev. D 54, (1996). 

P. Anninos, D. Hobill, E. Seidel, L. Smarr, and W.-M. 

Suen, Phys. Rev. Lett. 71, 2851 (1993). 

P. Anninos, D. Hobill, E. Seidel, L. Smarr, and W.-M. 

Suen, Phys. Rev. D 52, 2044 (1995). 

K. Camarda, Ph.D. thesis, University of Illinois at 

Urbana-Champaign, Urbana, Illinois, 1998. 

K. Ca marda and E. Seidel, Phys. Rev. D 57, 3204 

(1998), |gr-qc/9709075j 

G. Allen, K. Camarda, and E. Seidel, (1998), in prepa- 
ration. 

S. L. Shapiro and S. A. Teukolsky, in Dynamical Space- 
times and Numerical Relativity, edited by J. M. Cen- 
trella (Cambridge University Press, Cambridge, Eng- 
land, 1986), pp. 74-100. 

R. H. Price and J. Pullin, Phys. Rev. Lett. 72, 3297 

(1994) . 

P. Anninos, R. H. Price, J. Pullin, E. Seidel, and W.-M. 

Suen, Phys. Rev. D 52, 4462 (1995). 

A. Abrahams and R. Price, Phys. Rev. D 53, 1972 

(1996). 

J. Baker, A. Abrahams, P. Anninos, S. Brandt, R. Price, 
J. Pullin, and E. Seidel, Phys. Rev. D 55, 829 (1997). 
R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, 
Physical Review Letters 77, 4483 (1996). 
L. Rezzolla, A. M. Abrahams, T. W. Baumgarte, G. B. 
Cook, M. A. Scheel, S. L. Shapiro, and S. A. Teukolsky, 
Phys. Rev. D 57, 1084 (1998). 

M. Berger and J. Oliger, Journal of Computational 



[34 
[35; 

[36; 

[37 

[38; 

[39; 

[40; 

[41 

[42; 
[43; 



[44 
[45 

[46" 
[47 



[49" 
[50" 
[51 

[52; 

[53" 
[54 
[5& 
[56 

[57 



[5s; 
[59; 

[60; 

[61 
[62; 
[63; 
[64 



Physics 53, 484 (1984). 

L. A. Wild, Ph.D. thesis, University of Wales, 1996. 
P. Papadapoulos, E. Seide l, and L. Wild , Physical Re- 
view D (1998), submitted, ^r-qc/980206£ . 
E. Seidel and W.-M. Suen, Phys. Rev. Lett. 69, 1845 
(1992). 

P. Anninos, G. Daues, J. Masso, E. Seidel, and W.-M. 
Suen, Phys. Rev. D 51, 5562 (1995). 
M. A. Scheel, S. L. Shapiro, and S. A. Teukolsky, Phys. 
Rev. D 51, 4208 (1995). 

R. Marsa and M. Choptuik, Phys Rev D 54, 4929 
(1996). 

G. E. Daues, Ph.D. thesis, Washington University, St. 
Louis, Missouri, 1996. 

P. Anninos, J. Masso, E. Seidel, W.-M. Suen, and M. 

Tobias, Phys. Rev. D 56, 842 (1997). 

P. Anninos, J. Masso, E. Seidel, W.-M. Suen, and M. 

Tobias, Phys. Rev. D 54, 6544 (1996). 

J. Balakrishna, G. Daues, E. Seidel, W.-M. Suen, M. 

Tobias, and E. Wang, Class. Quant. Grav. 13, L135 

(1996). 

M. Alcubierre, Phys. Rev. D 55, 5981 (1997). 

M. Alcubierre an d J. Masso, P hys. Rev. D Rapid 



Comm., to appear, gr-qc/9709024 . 

B. Briigmann, (1997), |gr-qc/9708035 . 

R. Gomez, L. Lehner, R. Marsa, and J. Winicour, 



(1997), |gr-qc/9710138 . 

J. York, in Sources of Gravitational Radiation, edited 
by L. Smarr (Cambridge University Press, Cambridge, 
England, 1979). 

H. Friedrich, Class. Quant. Grav. 13, 1451 (1996). 

O. Reula, Living Reviews in Relativity 1, (1998). 

R. J. Leveque, Numerical Methods for Conservation 

Laws (Birkhauser Verlag, Basel, 1992). 

Y. Choquet-Bruhat and T. Ruggeri, Comm. Math. Phys 

89, 269 (1983). 

H. Friedrich, Comm. Math. Phys. 100, 1525 (1985). 
C. Bona and J. Masso, Phys. Rev. D 38, 2419 (1988). 
C. Bona and J. Masso, Phys. Rev. D 40, 1022 (1989). 
J. Masso, Ph.D. thesis, University of the Balearic Is- 
lands, 1992. 

C. Bona and J. Masso, in Approaches to Numerical Rel- 
ativity, edited by R. D'Inverno (Cambridge University 
Press, Cambridge, England, 1992). 

R. Gjertsen, J. Masso, M. Nardulli, E. Seidel, J. Shalf, 
and D. Weber, Forefronts 11, (1996), Cornell Theory 
Center. 

S. Kuo, M. Winslett, K. Seamons, Y. Chen, Y. Cho, 
and M. Subramaniam, in Proceedings of the SIAM Con- 
ference on Parallel Processing for Scientific Computing 
(SIAM, Minneapolis, MN, 1997). 

C. Bona, J. Masso, and J. Stela, Phys. Rev. D 51, 1639 
(1995). 

A. Arbona, C. Bona, J. Masso, and J. Stela, In prepa- 
ration (1998). 

S. Fritelli and O. Reula, Commun. Math. Phys. 166, 
221 (1994). 

Y. Choquet-Bruhat and J. York, C. R. Acad. Sc. Paris 
321, 1089 (1995). 

A. Abrahams, A. Anderson, Y. Choquet-Bruhat, and J. 



34 



York, Phys. Rev. Lett. 75, 3377 (1995). 

S. Fritelli and O. Reula, Phys. Rev. Lett. 75, 4667 

(1995). 

M. H. van Putten and D. Eardley, Phys. Rev. D 53, 
3056 (1996). 

A. Abrahams, A. Anderson, Y. Choquet-Bruhat, and J. 
York, Class.Quant.Grav A9 (1997). 
M. Scheel, T. Baumgarte, G. Cook, S. Shapiro, and S. 
Teukolsky, Phys. Rev. D 56, 6320 (1997). 
G. Cook and M. Scheel, private communication. 
T. Clune, J. Masso, M. Miller, and P. Walker, Techni- 
cal report, National Center for Supercomputing Appli- 
cations, (unpublished), in preparation. 
P. Walker, Technical report, National Center for Super- 
computing Applications, (unpublished), in preparation. 
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. 
Vetterling, Numerical Recipes (Cambridge University 
Press, Cambridge, England, 1986). 

A. Abrahams, A. Anderson, Y. Choquet-Bruhat, and J. 

York, C.R. Acad. Sci. Paris 835 (1996). 

The Message Passing Inte rface (MPI) standard, 



http:/ /www. ncs.anl.gov/mpi/ 



M. Parashar and J. C. Browne, Proceedings of the Inter- 
national Conference for High Performance Computing 
22 (1995). 

M. Parashar and J. C. Browne, Proceedings of the 29 th 
Annual Hawaii International Conference on System Sci- 
ences 604 (1996). 

P. Anninos, D. Hobill, E. Seidel, and W.-M. Suen, Pro- 
ceedings of the Sixth Canadian Conference on General 
Relativity and Relativistic Astrophysics (Fields Institute 
Comunications, 1996), in press. 

G. Allen, K. Camarda, and E. Seidel, in preparation. 
N. Bishop, R. Gomez, L. Lehner, and J. Winicour, Phys. 
Rev. D 52, (1997). 



P. Hiibner, ;r-qc/980406£ (1998) 



P. Hiibner, Phys. Rev. D 53, 701 (1996). 

J. Frauendiener, gr-qc/971250 and |gr-qc/9712052 
(1997). 

A. Abrahams, D. Bernstein, D. Hobill, E. Seidel, and L. 

Smarr, Phys. Rev. D 45, 3544 (1992). 

A. Arbona, C. Bona, J. Carot, L. Mas, J. Masso, and 

J. Stela, Phys. Rev. D 57, 2397 (1998). 

C. Bona, in Relativity and Scientific Computing, edited 

by F. Hehl (Springer- Verlag, Berlin, 1996). 

H. C. Yee, Computational Fluid Dynamics (von Karman 

Institute for Fluid Dynamics, NASA Ames Research 

Center, CA, 1989). 

P. Walker, Ph.D. thesis, University of Illinois at Urbana- 
Champaign, Urbana, Illinois, 1998, in preparation. 
M. Choptuik, Phys. Rev. D 44, 3124 (1991). 
L. Gunnarsen, H. Shinkai, and K. Maeda, Class. Quant. 
Grav. 12, 133 (1995). 

P. Anninos, K. Camarda, J. Masso, E. Seidel, W.-M. 
Suen, M. Tobias, and J. Towns, in The Seventh Mar- 
cel Grossmann Meeting: On Recent Developments in 
Theoretical and Experimental General Relativity, Grav- 
itation, and Relativistic Field Theories, edited by R. T. 
Jantzen, G. M. Keiser, and R. Ruffini (World Scientific, 
Singapore, 1996), pp. 644-647. 



[91 

[92; 
[93; 

[94" 
[95 



[96" 
[97 



[99; 

[100; 

[101 
[102; 



[103; 
[104 



U. Yurtsever, Phys. Rev. D 37, 2790 (1988). 

U. Yurtsever, Phys. Rev. D 38, 1731 (1988). 

M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 

(1995). 

S. Teukolsky, Phys. Rev. D 26, 745 (1982). 

C. Evans, in Dynamical Spacetimes and Numerical Rel- 
ativity, edited by J. Centrella (Cambridge University 
Press, Cambridge, England, 1986), pp. 3-39. 

K. Camarda and E. Seidel, in preparation. 

D. Bernstein, D. Hobill, E. Seidel, L. Smarr, and J. 
Towns, Phys. Rev. D 50, 5000 (1994). 

S. Brandt and E. Seidel, Phys. Rev. D 52, 856 (1995). 
S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 
(1997). 

C. W. Misner, K. S. Thorne, and J. A. Wheeler, Grav- 
itation (W. H. Freeman, San Francisco, 1973). 
B. Briigmann, private communication. 
S. Balay, W. Gropp, L. C. Mclnnes, and B. Smith, 
PETSc - The Portable, Extensible Toolkit for Scientific 



Compuation, 1998, http: //www 
D. Bernstein, Ph.D. thesis 



mcs.anl.gov/petsc/ 



University of Illinois 

Urbana-Champaign, 1993. 

P. Anninos, D. Bernstein, D. Hobill, E. Seidel, L. Smarr, 
and J. Towns, in Computational Astrophysics: Gas Dy- 
namics and Particle Methods, edited by W. Benz, J. 
Barnes, E. Muller, and M. Norman (Springer- Verlag, 
New York, 1997), in press. 
[105] C. Gundlach, Phys. Rev. D 57, 863 (1998). 



35 



